Commit e7646339 authored by s.zitz's avatar s.zitz
Browse files

Update of docstrings

parent 75132759
Pipeline #22916 passed with stage
in 5 minutes and 35 seconds
......@@ -145,6 +145,7 @@ notable changes.
- 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.
- Docstrings of all functions have been updated accordingly.
### Removed
......@@ -152,7 +153,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).
- `filmpressure!` in favor of `h∇p!`, see [forcing.jl](src/forcing.jl).
## [0.0.1] - 2020-08-18
......
......@@ -106,9 +106,9 @@ version = "0.10.9"
[[ColorVectorSpace]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "StatsBase"]
git-tree-sha1 = "bd0c0c81a39923bc03f9c3b61d89ad816e741002"
git-tree-sha1 = "2ae827d936fa9d8e00dd5166563499c07c5672c5"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.8.5"
version = "0.8.6"
[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"]
......@@ -124,9 +124,9 @@ version = "0.3.0"
[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "215f1c81cfd1c5416cd78740bff8ef59b24cd7c0"
git-tree-sha1 = "7c7f4cda0d58ec999189d70f5ee500348c4b4df1"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.15.0"
version = "3.16.0"
[[CompilerSupportLibraries_jll]]
deps = ["Libdl", "Pkg"]
......@@ -174,10 +174,10 @@ uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "0.21.7"
[[DataStructures]]
deps = ["InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "88d48e133e6d3dd68183309877eac74393daa7eb"
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "a88a67acbf3b61057371f315cadd027c8bce6d6d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.17.20"
version = "0.18.5"
[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
......@@ -490,10 +490,10 @@ uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.2"
[[JuMP]]
deps = ["Calculus", "DataStructures", "ForwardDiff", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "cbab42e2e912109d27046aa88f02a283a9abac7c"
deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "766014f271bd33b7f9d9bdc4847e214ee20ae84d"
uuid = "4076af6c-e467-56ae-b986-b466b2749572"
version = "0.21.3"
version = "0.21.4"
[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
......@@ -556,9 +556,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
[[MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "cd2049c055c7d192a235670d50faa375361624ba"
git-tree-sha1 = "a121678315668752aa5d26fef9ab564941c81e43"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.9.14"
version = "0.9.15"
[[MathProgBase]]
deps = ["LinearAlgebra", "SparseArrays"]
......@@ -611,9 +611,9 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.4"
[[OffsetArrays]]
git-tree-sha1 = "b8500f9d73999cfbab4add5136ec26894081581e"
git-tree-sha1 = "663d3402efa943c95f4736fa7b462e9dd97be1a9"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.1.3"
version = "1.2.0"
[[OpenSpecFun_jll]]
deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
......@@ -622,9 +622,9 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.3+3"
[[OrderedCollections]]
git-tree-sha1 = "293b70ac1780f9584c89268a6e2a560d938a7065"
git-tree-sha1 = "16c08bf5dba06609fe45e30860092d6fa41fde7b"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.3.0"
version = "1.3.1"
[[PaddedViews]]
deps = ["OffsetArrays"]
......
......@@ -13,7 +13,7 @@ include("runsimulation.jl")
include("measurements.jl")
# Export generators for structs
export default_moments, default_cumoments, default_velocity, default_cuvelocity, default_force, default_cuforce
export default_moments, default_cumoments, default_velocity, default_cuvelocity, default_force, default_cuforce, default_distribution
# Defined initial fluid states
export flatsurface, sinewaves, singledroplet
# As well as initial substrate states
......
......@@ -6,6 +6,7 @@ Computes the macroscopic moments `height` and `velocity` from the post collision
In case of using `BGKwithBuickGreatedStream` or `BGKwithGuoStream` add the sum of all forces to the moment calculation.
# Theory
The macroscopic moments, `density` (ρ) and `velocity` can be obtained as moments of the distribution function.
Density ρ reduces in the shallow water case to the `height` (h) of the film and is given by the zeroth moment.
On the other hand the velocity vector in 2d `(u,v)` is given by both the first and the second moment.
......@@ -23,10 +24,11 @@ To be efficient the height is computed simply using **Julia**s `sum!` function.
For the velocities I decided to do a little index shifting to use `sum` as well.
# Examples
```jldoctest
julia> using JuThinFilm
julia> dists = distributions{Float64}(para=params(lx=5,ly=5));
julia> dists = default_distribution(Float64, 5, 5);
julia> mom = default_moments(Float64, 5, 5); mom.height
5×5 Array{Float64,2}:
......@@ -48,7 +50,7 @@ julia> mom.height # height is now the sum of fout.
0.1 0.1 0.1 0.1 0.1
0.1 0.1 0.1 0.1 0.1
julia> dists = cudistributions{Float32}(para=params{Float32}(lx=5,ly=5)); dists.fout[:, :, 1] .= 0.1f0; # Works similar on the GPU.
julia> dists = JuThinFilm.cudistributions{Float32}(para=JuThinFilm.params{Float32}(lx=5,ly=5)); dists.fout[:, :, 1] .= 0.1f0; # Works similar on the GPU.
julia> mom = default_cumoments(Float32, 5, 5); mom.height
5×5 CUDA.CuArray{Float32,2}:
......
......@@ -7,6 +7,7 @@ This is a simple BKG collision kernel with additional streaming done by `circshi
To include the forcing the weighting factor method (WFM) is used.
# Theory
The idea behind this collision operator is quite simple.
Within the time `τ` a system close to equilibrium should relax towards its equilibrium.
......@@ -22,6 +23,7 @@ Such the computation of ``f_{\\alpha}(\\mathbf{x}+\\mathbf{c}_{\\alpha}\\Delta t
This is straightforwardly taken care of by Julias function `circshift`, which has the additional feature of computing `ftemp` for the next time step.
# Arguments
- `fout::Array{T,3}`: output distribution function after the collision process.
- `ftemp::Array{T,3}`: distribution function from the previous time step.
- `feq::Array{T,3}`: equilibrium distribution function.
......@@ -32,10 +34,11 @@ This is straightforwardly taken care of by Julias function `circshift`, which ha
# Examples
```jldoctest
julia> using JuThinFilm, Test
julia> dists = distributions{Float64}(para=params(lx=5,ly=5)); # allocate some distributions
julia> dists = default_distribution(Float64, 5, 5); # allocate some distributions
julia> dists.feq[:,:,1] = ones(5,5)
5×5 Array{Float64,2}:
......@@ -54,26 +57,27 @@ julia> dists.fout[:,:,1]
0.0 0.0 0.0 0.0 0.0
julia> f = default_force(Float64, 5, 5)
force{Float64}
para: params{Float64}
JuThinFilm.force{Float64}
para: JuThinFilm.params{Float64}
x: Array{Float64}((5, 5)) [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]
y: Array{Float64}((5, 5)) [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> col = collisionparams{Float64}(); # default τ=1
julia> col = JuThinFilm.collisionparams(); # default τ=1
julia> JuThinFilm.BGKandStream!(dists.fout, dists.ftemp, dists.feq, f.x, f.y, col)
julia> BGKandStream!(dists.fout, dists.ftemp, dists.feq, f.x, f.y, col)
julia> @test all(dists.fout == dists.feq) # With τ=1 and no forcing fout is similar to feq!
Test Passed
julia> dists.fout .= 0.0
julia> dists.fout .= 0.0;
julia> JuThinFilm.BGKandStream!(dists, f, col)
julia> BGKandStream!(dists, f, col)
julia> @test all(dists.fout == dists.feq) # With τ=1 and no forcing fout is similar to feq!
Test Passed
```
# References
The omnipotent reference on the issue of the lattice Boltzmann method is the book by Timm Krüger
- [The lattice Boltzmann Method](https://www.springer.com/gp/book/9783319446479)
......@@ -86,10 +90,10 @@ Specially for the shallow water lattice Boltzmann
The last paper is our own, it summarizes nicely how we solve the problems at hand.
See also: [`circshift`](@ref), [`runsimulation`](@ref), [`BGKwithGUOcorrectionandStream!`](@ref)
See also: [`circshift`](@ref), [`runsimulation`](@ref), [`BGKwithBuickGreatedStream!`](@ref)
"""
function BGKandStream!(fout, feq, ftemp, Fx, Fy, col::collisionparams; kind="simple")
function BGKandStream!(fout, feq, ftemp, Fx, Fy, col::JuThinFilm.collisionparams; kind="simple")
# The lattice speeds
cx = [0 1 0 -1 0 1 -1 -1 1]
cy = [0 0 1 0 -1 1 1 -1 -1]
......@@ -113,7 +117,7 @@ function BGKandStream!(fout, feq, ftemp, Fx, Fy, col::collisionparams; kind="sim
return nothing
end
# Dispatch with predefined mutable structs!
function BGKandStream!(dists::Distribution, F::Force, col::collisionparams; kind="simple")
function BGKandStream!(dists::JuThinFilm.Distribution, F::JuThinFilm.Force, col::JuThinFilm.collisionparams; kind="simple")
# The lattice speeds
cx = [0 1 0 -1 0 1 -1 -1 1]
cy = [0 0 1 0 -1 1 1 -1 -1]
......@@ -145,6 +149,7 @@ Pretty much the same as `BGKandStream!` but with a different force correction.
If used half of the forcing has to be added to the moment calculation.
# Theory
The only difference is in the definition of the force correction term.
``\\mathcal{S}_{\\alpha} = (1 - \\frac{1}{2\\tau})\\frac{3 w_{\\alpha}\\mathbf{c}_{\\alpha}}{c^2}\\cdot\\mathbf{F}.``
......@@ -154,6 +159,7 @@ Comparing those two, the `WFM` picks up a factor of 0.5 for the diagonal terms.
Indeed for the sliding droplet the `WFM` approach yields the correct droplet shape while the `simple` one comes short.
# Arguments
- `fout::Array{T,3}`: output distribution function after the collision process.
- `ftemp::Array{T,3}`: distribution function from the previous time step.
- `feq::Array{T,3}`: equilibrium distribution function.
......@@ -163,10 +169,11 @@ Indeed for the sliding droplet the `WFM` approach yields the correct droplet sha
- `kind::String="simple"`: force correction, either "simple" or "WFM".
# Examples
```jldoctest
julia> using JuThinFilm, Test
julia> dists = distributions{Float64}(para=params(lx=5,ly=5)); # allocate some distributions
julia> dists = default_distribution(Float64, 5, 5); # allocate some distributions
julia> dists.feq[:,:,1] = ones(5,5)
5×5 Array{Float64,2}:
......@@ -185,26 +192,27 @@ julia> dists.fout[:,:,1]
0.0 0.0 0.0 0.0 0.0
julia> f = default_force(Float64, 5, 5)
force{Float64}
para: params{Float64}
JuThinFilm.force{Float64}
para: JuThinFilm.params{Float64}
x: Array{Float64}((5, 5)) [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]
y: Array{Float64}((5, 5)) [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> col = collisionparams{Float64}(); # default τ=1
julia> col = JuThinFilm.collisionparams{Float64}(); # default τ=1
julia> JuThinFilm.BGKwithBuickGreatedStream!(dists.fout, dists.ftemp, dists.feq, f.x, f.y, col)
julia> BGKwithBuickGreatedStream!(dists.fout, dists.ftemp, dists.feq, f.x, f.y, col)
julia> @test all(dists.fout == dists.feq) # With τ=1 and no forcing fout is similar to feq!
Test Passed
julia> dists.fout .= 0.0
julia> dists.fout .= 0.0;
julia> JuThinFilm.BGKwithBuickGreatedStream!(dists, f, col)
julia> BGKwithBuickGreatedStream!(dists, f, col)
julia> @test all(dists.fout == dists.feq) # With τ=1 and no forcing fout is similar to feq!
Test Passed
```
# References
- [Body Force Term Comparison in Lattice Boltzmann Method for Shallow Water Equations](https://aip.scitation.org/doi/pdf/10.1063/1.4825929)
- [Gravity in a lattice Boltzmann model](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.61.5307)
......@@ -214,7 +222,7 @@ Therefore I assume we are fine.
See also: [`circshift`](@ref), [`runsimulation`](@ref), [`BGKandStream!`](@ref), [`BGKwithGuoStream`](@ref)
"""
function BGKwithBuickGreatedStream!(fout, feq, ftemp, Fx, Fy, col::collisionparams; kind="simple")
function BGKwithBuickGreatedStream!(fout, feq, ftemp, Fx, Fy, col::JuThinFilm.collisionparams; kind="simple")
# The lattice speeds and weights
cx = [0 1 0 -1 0 1 -1 -1 1]
cy = [0 0 1 0 -1 1 1 -1 -1]
......@@ -239,7 +247,7 @@ function BGKwithBuickGreatedStream!(fout, feq, ftemp, Fx, Fy, col::collisionpara
return nothing
end
function BGKwithBuickGreatedStream!(dists::Distribution, F::Force, col::collisionparams; kind="simple")
function BGKwithBuickGreatedStream!(dists::JuThinFilm.Distribution, F::JuThinFilm.Force, col::JuThinFilm.collisionparams; kind="simple")
# The lattice speeds and weights
cx = [0 1 0 -1 0 1 -1 -1 1]
cy = [0 0 1 0 -1 1 1 -1 -1]
......@@ -278,25 +286,27 @@ If used half of the forcing has to be added to the moment calculation.
Which is taken care of if a force `F` is given as argument to `macroscopicmoments!`.
# Theory
The only difference is in the definition of the force correction term.
``\\mathcal{S}_{\\alpha} = (1 - \\frac{1}{2\\tau})w_{\\alpha}\\Bigg[3\\frac{\\mathbf{c}_{\\alpha}-\\mathbf{u}}{c^2}+9\\frac{\\mathbf{c}_{\\alpha}\\mathbf{u}}{c^2}\\mathbf{c}_{\\alpha}\\Bigg]\\cdot\\mathbf{F}.``
# Arguments
- `fout::Array{T,3}`: output distribution function after the collision process.
- `ftemp::Array{T,3}`: distribution function from the previous time step.
- `feq::Array{T,3}`: equilibrium distribution function.
- `Fx::Array{T,2}`: x-component of forces applied to the fluid.
- `Fy::Array{T,2}`: y-component of forces applied to the fluid.
- `col::collisionparams`: contains `τ` and `ω`.
- `kind::String="simple"`: force correction, either "simple" or "WFM".
- `kind::String`: force correction, either "simple" or "WFM" default is set to "simple".
# Examples
```jldoctest
julia> using JuThinFilm, Test
julia> dists = distributions{Float64}(para=params(lx=5,ly=5)); # allocate some distributions
julia> dists = default_distribution(Float64, 5, 5); # allocate some distributions
julia> dists.feq[:,:,1] = ones(5,5)
5×5 Array{Float64,2}:
......@@ -315,21 +325,23 @@ julia> dists.fout[:,:,1]
0.0 0.0 0.0 0.0 0.0
julia> f = default_force(Float64, 5, 5)
force{Float64}
para: params{Float64}
JuThinFilm.force{Float64}
para: JuThinFilm.params{Float64}
x: Array{Float64}((5, 5)) [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]
y: Array{Float64}((5, 5)) [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> col = collisionparams{Float64}(); # default τ=1
julia> col = JuThinFilm.collisionparams{Float64}(); # default τ=1
julia> mom = default_moments(Float64,5,5); mom.velocity.x .= 0.1; mom.velocity.y .= -0.1;
julia> JuThinFilm.BGKwithBuickGreatedStream!(dists.fout, dists.ftemp, dists.feq, f.x, f.y, col)
julia> BGKwithGuoStream!(dists.fout, dists.ftemp, dists.feq, mom.velocity.x, mom.velocity.y, f.x, f.y, col)
julia> @test all(dists.fout == dists.feq) # With τ=1 and no forcing fout is similar to feq!
Test Passed
julia> dists.fout .= 0.0
julia> dists.fout .= 0.0;
julia> JuThinFilm.BGKwithBuickGreatedStream!(dists, f, col)
julia> BGKwithGuoStream!(dists, mom, f, col)
julia> @test all(dists.fout == dists.feq) # With τ=1 and no forcing fout is similar to feq!
Test Passed
......@@ -340,7 +352,7 @@ Test Passed
See also: [`circshift`](@ref), [`runsimulation`](@ref), [`BGKandStream!`](@ref), [`BGKwithBuickGreatedStream`](@ref), [`macroscopicmoments!`](@ref)
"""
function BGKwithGuoStream!(fout, feq, ftemp, velx, vely, Fx, Fy, col::collisionparams; kind="simple")
function BGKwithGuoStream!(fout, feq, ftemp, velx, vely, Fx, Fy, col::JuThinFilm.collisionparams; kind="simple")
# The lattice speeds and weights
cx = [0 1 0 -1 0 1 -1 -1 1]
cy = [0 0 1 0 -1 1 1 -1 -1]
......
......@@ -4,6 +4,7 @@
Computes the equilibrium distribution `dists.feq` of a given macroscopic state `mom`.
# Theory
One of the reasons why the lattice Boltzmann method became so widely spread is due to its simple assumptions.
Among them is that any system simulated with lattice Boltzmann is a system close to equilibrium.
For an ideal gas one knows that the momentum distribution is given by a maxwellian distribution with temperature T,
......@@ -28,20 +29,21 @@ as zeroth population and
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.
- `lat::lattice`: struct that stores the lattice velocities `ci` and weights `wi`.
- `parameters::param`: struct that stores lattice size and gravitational acceleration.
- `dists::JuThinFilm.distributions`: container for the three distribution functions.
- `mom::JuThinFilm.moments`: struct that contains a macroscopic state, e.g. height and velocity.
- `lat::JuThinFilm.lattice`: struct that stores the lattice velocities `ci` and weights `wi`.
- `parameters::JuThinFilm.param`: struct that stores lattice size and gravitational acceleration.
# Examples
```jldoctest
julia> using JuThinFilm
julia> dists = JuThinFilm.distributions{Float64}(para = params(lx=5, ly=5));
julia> dists = default_distribution(Float64, 5, 5);
julia> mom = JuThinFilm.default_moments(Float64, 5, 5);
julia> mom = default_moments(Float64, 5, 5);
julia> lat = lattice{Float64}(); par = params{Float64}(lx=5, ly=5);
julia> lat = JuThinFilm.lattice{Float64}(); par = JuThinFilm.params{Float64}(lx=5, ly=5);
julia> equilibrium!(dists, mom, lat, par)
......@@ -53,9 +55,9 @@ julia> dists.feq[:, :, 1] # Now it should be filled with mom.height
1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0
julia> mom.velocity.x = fill(0.1, (5,5)); # With some small velocity
julia> mom.velocity.x .= 0.1; # With some small velocity
julia> equilibrium!(dists, mom, lat, par, Float64)
julia> equilibrium!(dists, mom, lat, par)
julia> dists.feq[:, :, 1]
5×5 Array{Float64,2}:
......@@ -121,13 +123,13 @@ Computes the square of the velocity vector `mom.velocity`, where `mom` is a mome
```jldoctest
julia> using JuThinFilm
julia> par = params{Float64}(lx=5, ly=5); vel = velocity{Float64}(para = par);
julia> par = JuThinFilm.params{Float64}(lx=5, ly=5); vel = JuThinFilm.velocity{Float64}(para = par);
julia> mom = moments{Float64}(para = par, velocity = vel)
moments{Float64}
para: params{Float64}
julia> mom = default_moments(Float64,5,5)
JuThinFilm.moments{Float64}
para: JuThinFilm.params{Float64}
height: Array{Float64}((5, 5)) [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]
velocity: velocity{Float64}
velocity: JuThinFilm.velocity{Float64}
pressure: Array{Float64}((5, 5)) [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> JuThinFilm.velsquare(mom) # Because vel is just 0s
......@@ -138,9 +140,9 @@ julia> JuThinFilm.velsquare(mom) # Because vel is just 0s
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
julia> vel2 = velocity{Float64}(para = par, x=fill(2.0, 5, 5), y=fill(2.0, 5, 5));
julia> vel2 = default_velocity(Float64, 5, 5); vel2.x .= 2; vel2.y .= 2;
julia> mom2 = moments{Float64}(para = par, velocity = vel2);
julia> mom2 = JuThinFilm.moments{Float64}(para = par, velocity = vel2);
julia> JuThinFilm.velsquare(mom2) # Now it should be 2*2 + 2*2 = 8
5×5 Array{Float64,2}:
......
"""
slippage!(fx, fy, h, vx, vy, col, T)
slippage!(fx, fy, h, vx, vy, col)
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.
......@@ -23,23 +24,24 @@ The force should be strong if `` h \\rightarrow 0 ``, therefore the full term is
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.
- `fy::Array{T,2}`: y-component of the computed force.
- `h::Array{T,2}`: height field at every grid point.
- `vx::Array{T,2}`: x-component of the velocity.
- `vy::Array{T,2}`: y-component of the velocity.
- `col::collisionparams`: struct that contains slip lenght `δ` and viscosity `μ`.
- `T:: <:AbstractFloat`: numerical output type of result, e.g. `Float64`.
- `col::JuThinFilm.collisionparams`: struct that contains slip lenght `δ` and viscosity `μ`.
# Examples
```jldoctest
julia> using JuThinFilm
julia> f = force{Float64}(para=params(lx=5,ly=5)); mom = default_moments(Float64, 5, 5);
julia> f = default_force(Float64, 5, 5); mom = default_moments(Float64, 5, 5);
julia> col = collisionparams{Float64}();
julia> col = JuThinFilm.collisionparams();
julia> slippage!(f.x, f.y, mom.height, mom.velocity.x, mom.velocity.y, col, Float64) # No velocity and no unit height f=0.
julia> slippage!(f.x, f.y, mom.height, mom.velocity.x, mom.velocity.y, col) # No velocity and no unit height f=0.
julia> f.x
5×5 Array{Float64,2}:
......@@ -49,7 +51,7 @@ julia> f.x
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
julia> mom.velocity.x .= 0.1; slippage!(f.x, f.y, mom.height, mom.velocity.x, mom.velocity.y, col, Float64)
julia> mom.velocity.x .= 0.1; slippage!(f.x, f.y, mom.height, mom.velocity.x, mom.velocity.y, col)
julia> f.x
5×5 Array{Float64,2}:
......@@ -59,11 +61,13 @@ julia> f.x
0.00909091 0.00909091 0.00909091 0.00909091 0.00909091
0.00909091 0.00909091 0.00909091 0.00909091 0.00909091
julia> slippage!(f, mom, col, Float64) # Can also be called with structs
julia> f.x .= 0.0; f.y .= 0.0;
julia> slippage!(f, mom, col) # Can also be called with structs
julia> cumom = default_cumoments(Float32, 5, 5); cuf = cuforce{Float32}(para = params{Float32}(lx=5,ly=5));
julia> cumom = default_cumoments(Float32, 5, 5); cuf = default_cuforce(Float32, 5, 5); col32 = JuThinFilm.collisionparams{Float32}();
julia> cumom.velocity.x .= 0.1f0; slippage!(cuf, cumom, col, Float32) # And on the GPU
julia> cumom.velocity.x .= 0.1f0; slippage!(cuf, cumom, col32) # And on the GPU
julia> cuf.x
5×5 CUDA.CuArray{Float32,2}:
......@@ -83,7 +87,7 @@ There is plenty of literature dealing with it, Fetzer and Jacobs are definitely
The second reference is the paper about the method.
See also: [`collidestreamBGK_periodic!`](@ref)
See also: [`BGKandStream!`](@ref), [`runsimulation`](@ref)
"""
function slippage!(Fx, Fy, height, velocityx, velocityy, col::collisionparams)
......@@ -144,15 +148,13 @@ julia> using JuThinFilm, Random, Statistics, Test
julia> lx = 100; ly = 100;
julia> f = force{Float64}(para=params(lx=lx,ly=ly));
julia> r = force{Float64}(para=params(lx=lx,ly=ly));
julia> f = default_force(Float64, lx, ly);
julia> mom = default_moments(Float64, lx, ly);
julia> col = collisionparams{Float64}(kᵦT = 11/6, μ=2); # Generates a nice variance
julia> col = JuThinFilm.collisionparams{Float64}(kᵦT = 11/6, μ=2); # Generates a nice variance
julia> thermalflutuations!(f, mom, r, col, Float64)
julia> thermalfluctuations!(f, mom, col)
julia> @test sum(f.x)/(lx*ly) ≈ 0.0 atol = 1e-1 # check the mean
Test Passed
......@@ -169,9 +171,9 @@ A few good articles to get a general understanding
- [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)
See also: [`BGKandStream!`](@ref)
"""
function thermalfluctuations!(Fx, Fy, height, col::collisionparams)
function thermalfluctuations!(Fx, Fy, height, col::JuThinFilm.collisionparams)
# Some random number generation, have to switch between GPU and CPU
if isa(Fx, CuArray)
Fx .= CUDA.randn(Float64, size(height))
......@@ -193,7 +195,7 @@ function thermalfluctuations!(Fx, Fy, height, col::collisionparams)
return nothing
end
function thermalfluctuations!(F::Force, mom::Moment, col::collisionparams)
function thermalfluctuations!(F::JuThinFilm.Force, mom::JuThinFilm.Moment, col::JuThinFilm.collisionparams)
# Some random number generation, have to switch between GPU and CPU
randn!(F.x)
randn!(F.y)
......@@ -218,6 +220,7 @@ 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.
For the next iteration of this program add ghost cells to compute the boundary conditions more easily.
# Theory
......@@ -228,18 +231,18 @@ Starting with the thin film equation
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.
For details take a look at the [`pressure.jl`](src/pressure.jl) file.
# Examples
```jldoctest
julia> using JuThinFilm, CUDA
julia> using JuThinFilm, Images
julia> mom = default_moments(Float64,5,5); f = force{Float64}(para = params(lx=5,ly=5));
julia> mom = default_moments(Float64,5,5); f = default_force(Float64, 5, 5);
julia> mom.pressure .= 1.0; mom.height .= 0.5; # Constant pressure yields no gradient!
julia> h∇p!(f, mom)
julia> h∇p!(f, mom, padarray(mom.pressure, Pad(:circular, 1, 1)))
julia> f.x
5×5 Array{Float64,2}:
......@@ -257,7 +260,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> h∇p!(f, mom)
julia> h∇p!(f, mom, padarray(mom.pressure, Pad(:circular, 1, 1)))
julia> f.x
5×5 Array{Float64,2}:
......@@ -265,7 +268,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, lx, ly -0.75 -0.75 -0.75 -0.75
-0.75 -0.75 -0.75 -0.75 -0.75
julia> f.y
5×5 Array{Float64,2}:
......@@ -283,7 +286,7 @@ For the gradient of a *D2Q9* lattice Boltzmann field I suggested to use the pape
- [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)
See also: [`BGKandStream!`](@ref)
"""
function h∇p!(fx, fy, height, pressure)
......@@ -321,7 +324,7 @@ function h∇p!(f::Force, mom::Moment, pressure)
end
"""
inclinedplane!(F::Force, mom::Moment, t::Int, mag, x::Bool, y::Bool; tonset=1000, tsmooth=10000)
inclinedplane!(F::JuThinFilm.Force, mom::JuThinFilm.Moment, t::Int, mag, x::Bool, y::Bool; tonset=1000, tsmooth=10000)
Forcing as if the film was placed on an inclined plane, thus the heigher `height` the stronger the force.
......@@ -335,10 +338,10 @@ 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.
- `F::JuThinFilm.Force`: output forcing due to the chosen forcing strength `mag` and the height profile `mom.height`.
- `mom::JuThinFilm.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.
- `mag::AbstractFloat`: forcing strength, be aware to use small enough values e.g. `mag` ≪ 1.0!
- `mag<:AbstractFloat`: forcing strength, be aware to use small enough values e.g. `mag` ≪ 1.0!
- `x::Bool`: `true` if force is applied in x-direction.
- `y::Bool`: `true` if force is applied in y-direction.
- `tonset::Int`: **optional** time step when the force is starting to ramp up.
......@@ -352,13 +355,13 @@ julia> Fbody = default_force(Float64, 5, 5); mom = default_moments(Float64,5,5);
julia> inclinedplane!(Fbody, mom, 10, 1, true, false)