This should do for now

src/runsimulation.jl

-"""
-    runrealsimulation2(T, lx, ly, maxtime; kind="Droplet", θ=1/9, δ=4.0, h₀=25, θ₀=1/3)
+# """
+#     runrealsimulation2(T, lx, ly, maxtime; kind="Droplet", θ=1/9, δ=4.0, h₀=25, θ₀=1/3)
 
-Real simulation to test if a droplet is relaxing towards it's equilibrium contact angle.
-"""
-function runrealsimulation2(T, lx, ly, maxtime; kind="Droplet", θ=1/9, δ=4.0, h₀=25, θ₀=1/3)
-    center = (Int(lx/2), Int(ly/2))
-    if kind == "Droplet"
-        # mom = sinewaves(T, lx, ly)
-        mom = singledroplet(T, lx, ly, h₀, θ₀, center)
-    end
-    wetted = []
+# Real simulation to test if a droplet is relaxing towards it's equilibrium contact angle.
+# """
+# function runrealsimulation2(T, lx, ly, maxtime; kind="Droplet", θ=1/9, δ=4.0, h₀=25, θ₀=1/3)
+#     center = (Int(lx/2), Int(ly/2))
+#     if kind == "Droplet"
+#         # mom = sinewaves(T, lx, ly)
+#         mom = singledroplet(T, lx, ly, h₀, θ₀, center)
+#     end
+#     wetted = []
     
-    # Some unuseful IO
-    println("Hey mate, let's do some funky simulation!")
-    println("You choose to run on a simulation of a $kind")
-    # Some parameter defining at the start
-    parameters = params{T}(lx=lx, ly=ly, mt=maxtime)
-    collisionparameters = collisionparams{T}(δ=δ, τ=1.0)
-    pressureparameters = pressurestats{T}(θ=θ)
-    dists = distributions{T}(para=parameters)
-    lat = lattice{T}()
-    Fpressure = default_force(T, lx, ly)
-    Fsum = default_force(T, lx, ly)
-    Fslip = default_force(T, lx, ly)
-    Fbody = default_force(T, lx, ly)
-    JuThinFilm.measurecoxvoinov!(wetted, mom, collisionparameters, pressureparameters)
-    # Get the first equilibrium here
-    equilibrium!(dists, mom, lat, parameters, T)
-    copy!(dists.fout,dists.feq)
-    # Now move to the time loop
-    for t in 1:parameters.mt
-        if t % 100 == 0
-            mass = sum(mom.height)
-            println("time step $t with mass $mass")
-        end
-        copy!(dists.ftemp,dists.fout)
-        # Compute forces
-        kernel_pressure_2!(mom.pressure, mom.height, pressureparameters, T)
-        filmpressure!(Fpressure, mom, T)
-        slippage!(Fslip, mom, collisionparameters)
-        # constantforce!(Fbody, mom, t, 1e-4, true, false, T)
-        sumforces!(Fsum, Fpressure, Fslip, Fbody)
+#     # Some unuseful IO
+#     println("Hey mate, let's do some funky simulation!")
+#     println("You choose to run on a simulation of a $kind")
+#     # Some parameter defining at the start
+#     parameters = params{T}(lx=lx, ly=ly, mt=maxtime)
+#     collisionparameters = collisionparams{T}(δ=δ, τ=1.0)
+#     pressureparameters = pressurestats{T}(θ=θ)
+#     dists = distributions{T}(para=parameters)
+#     lat = lattice{T}()
+#     Fpressure = default_force(T, lx, ly)
+#     Fsum = default_force(T, lx, ly)
+#     Fslip = default_force(T, lx, ly)
+#     Fbody = default_force(T, lx, ly)
+#     JuThinFilm.measurecoxvoinov!(wetted, mom, collisionparameters, pressureparameters)
+#     # Get the first equilibrium here
+#     equilibrium!(dists, mom, lat, parameters, T)
+#     copy!(dists.fout,dists.feq)
+#     # Now move to the time loop
+#     for t in 1:parameters.mt
+#         if t % 100 == 0
+#             mass = sum(mom.height)
+#             println("time step $t with mass $mass")
+#         end
+#         copy!(dists.ftemp,dists.fout)
+#         # Compute forces
+#         kernel_pressure_2!(mom.pressure, mom.height, pressureparameters, T)
+#         filmpressure!(Fpressure, mom, T)
+#         slippage!(Fslip, mom, collisionparameters)
+#         # constantforce!(Fbody, mom, t, 1e-4, true, false, T)
+#         sumforces!(Fsum, Fpressure, Fslip, Fbody)
         
-        # Get the equilibrium and then performe a collision.
-        equilibrium!(dists, mom, lat, parameters)
-        JuThinFilm.collidestreamBGK_periodicGUO!(dists, Fsum, lat, collisionparameters, T)
-        JuThinFilm.macroscopicmoments!(mom, dists, Fsum)
-        JuThinFilm.measurecoxvoinov!(wetted, mom, collisionparameters, pressureparameters)
-    end
-    return mom, wetted
-end
-# That works seemingly well :)
-function slidingdroplet(T, lx, ly, maxtime; θ₀=1/3, θₑ=1/6, δ₀=2.5, R₀=40, fbody=8e-5)
-    center = (Int(lx/2), Int(ly/2))
-    mom = singledroplet(T, lx, ly, R₀, θ₀, center)
+#         # Get the equilibrium and then performe a collision.
+#         equilibrium!(dists, mom, lat, parameters)
+#         JuThinFilm.collidestreamBGK_periodicGUO!(dists, Fsum, lat, collisionparameters, T)
+#         JuThinFilm.macroscopicmoments!(mom, dists, Fsum)
+#         JuThinFilm.measurecoxvoinov!(wetted, mom, collisionparameters, pressureparameters)
+#     end
+#     return mom, wetted
+# end
+# # That works seemingly well :)
+# function slidingdroplet(T, lx, ly, maxtime; θ₀=1/3, θₑ=1/6, δ₀=2.5, R₀=40, fbody=8e-5)
+#     center = (Int(lx/2), Int(ly/2))
+#     mom = singledroplet(T, lx, ly, R₀, θ₀, center)
 
-    # Some parameter defining at the start
-    parameters = params{T}(lx=lx, ly=ly, mt=maxtime)
-    collisionparameters = collisionparams{T}(δ=δ₀, τ=1.0)
-    pressureparameters = pressurestats{T}(θ=θₑ, γ=0.007, hstar=0.1)
-    dists = distributions{T}(para=parameters)
-    lat = lattice{T}()
+#     # Some parameter defining at the start
+#     parameters = params{T}(lx=lx, ly=ly, mt=maxtime)
+#     collisionparameters = collisionparams{T}(δ=δ₀, τ=1.0)
+#     pressureparameters = pressurestats{T}(θ=θₑ, γ=0.007, hstar=0.1)
+#     dists = distributions{T}(para=parameters)
+#     lat = lattice{T}()
 
-    # Initialize Forcings
-    Fpressure = default_force(T, lx, ly)
-    Fsum = default_force(T, lx, ly)
-    Fslip = default_force(T, lx, ly)
-    Fbody = default_force(T, lx, ly)
-    hperiodic = zeros(lx+2, ly+2)
-    presperiodic = zeros(lx+2, ly+2)
+#     # Initialize Forcings
+#     Fpressure = default_force(T, lx, ly)
+#     Fsum = default_force(T, lx, ly)
+#     Fslip = default_force(T, lx, ly)
+#     Fbody = default_force(T, lx, ly)
+#     hperiodic = zeros(lx+2, ly+2)
+#     presperiodic = zeros(lx+2, ly+2)
 
-    # Some measuring lists
-    maxu = []
-    maxv = []
-    mass = []
-    # Get the first equilibrium here
-    JuThinFilm.equilibrium!(dists, mom, lat, parameters)
-    copy!(dists.fout,dists.feq)
-    for t in 1:parameters.mt
-        # push!(maxu, maximum(mom.velocity.x))
-        # push!(maxv, maximum(mom.velocity.y))
-        push!(mass, sum(mom.height))
-        if t%500  == 0 && t > 0
-            massout = round(mass[t], digits=3)
-            println("time step $t and mass $(massout)") 
-        end
-        # Compute forces
-        hperiodic = padarray(mom.height, Pad(:circular,1,1,))
-        JuThinFilm.capillarypressure!(mom.pressure, hperiodic, pressureparameters)
-        presperiodic = padarray(mom.pressure, Pad(:circular,1,1))
-        JuThinFilm.h∇p!(Fpressure.x, Fpressure.y, mom.height, presperiodic, lx, ly)
-        slippage!(Fslip, mom, collisionparameters)
-        # constantforce!(Fbody, mom, t, fbody, true, false, tonset=500, tsmooth=1000)
-        sumforces!(Fsum, Fpressure, Fslip, Fbody)
+#     # Some measuring lists
+#     maxu = []
+#     maxv = []
+#     mass = []
+#     # Get the first equilibrium here
+#     JuThinFilm.equilibrium!(dists, mom, lat, parameters)
+#     copy!(dists.fout,dists.feq)
+#     for t in 1:parameters.mt
+#         # push!(maxu, maximum(mom.velocity.x))
+#         # push!(maxv, maximum(mom.velocity.y))
+#         push!(mass, sum(mom.height))
+#         if t%500  == 0 && t > 0
+#             massout = round(mass[t], digits=3)
+#             println("time step $t and mass $(massout)") 
+#         end
+#         # Compute forces
+#         hperiodic = padarray(mom.height, Pad(:circular,1,1,))
+#         JuThinFilm.capillarypressure!(mom.pressure, hperiodic, pressureparameters)
+#         presperiodic = padarray(mom.pressure, Pad(:circular,1,1))
+#         JuThinFilm.h∇p!(Fpressure.x, Fpressure.y, mom.height, presperiodic, lx, ly)
+#         slippage!(Fslip, mom, collisionparameters)
+#         # constantforce!(Fbody, mom, t, fbody, true, false, tonset=500, tsmooth=1000)
+#         sumforces!(Fsum, Fpressure, Fslip, Fbody)
         
-        # Get the equilibrium and then performe a collision.
-        JuThinFilm.equilibrium!(dists, mom, lat, parameters)
-        JuThinFilm.BGKandStream!(dists.fout, dists.feq, dists.ftemp, Fsum.x, Fsum.y, collisionparameters, kind="WFM")
-        JuThinFilm.macroscopicmoments!(mom, dists)
-    end
-    # gif(anim, "animations/sliding_drop.gif", fps = 80)
-    return mom
-end
+#         # Get the equilibrium and then performe a collision.
+#         JuThinFilm.equilibrium!(dists, mom, lat, parameters)
+#         JuThinFilm.BGKandStream!(dists.fout, dists.feq, dists.ftemp, Fsum.x, Fsum.y, collisionparameters, kind="WFM")
+#         JuThinFilm.macroscopicmoments!(mom, dists)
+#     end
+#     # gif(anim, "animations/sliding_drop.gif", fps = 80)
+#     return mom
+# end

test/structs.jl

@@ -211,6 +211,23 @@
         @test size(mom.pressure) == (4,4)    
     end
 
+    @testset "default_distribution" begin
+        T = Float64
+        n, m  = 4, 4
+        dists = default_distribution(T, n, m)
+        @test isa(dists, JuThinFilm.distributions)
+        @test isa(dists.fout, Array{T, 3})
+        @test size(dists.fout) == (4,4,9)
+        @test all(dists.fout .== 0.0)
+        @test isa(dists.ftemp, Array{T, 3})
+        @test size(dists.ftemp) == (4,4,9)
+        @test all(dists.ftemp .== 0.0)   
+        @test isa(dists.feq, Array{T, 3})
+        @test size(dists.feq) == (4,4,9)
+        @test all(dists.feq .== 0.0)
+
+    end
+
     @testset "default_cumoments" begin
         T = Float32
         n, m  = 4, 4