WIP: GPU laplacian and gradient...

src/JuThinFilm.jl

 module JuThinFilm
 
-using CUDA, Revise, Parameters, Random, Images, LazySets, LinearAlgebra, ModularIndices
+using CUDA, Revise, Parameters, Random, Images, LazySets, LinearAlgebra, ModularIndices, Strided
 
 include("structs.jl")
 include("initialstates.jl")

src/pressure.jl

@@ -248,6 +248,29 @@ 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
+
+function capillarypressure2!(pressure::Array{Float64,2}, height, stats::pressurestats)
+    # The default compute size    
+    wid, len = size(pressure)
+    hᵪ = 0.05
+    h = @view height[2:len+1, 2:wid+1]
+    hip = @view height[3:len+2, 2:wid+1]
+    him = @view height[1:len, 2:wid+1]
+    hjp = @view height[2:len+1, 3:wid+2]
+    hjm = @view height[2:len+1, 1:wid]
+    hipjp = @view height[3:len+2, 3:wid+2]
+    himjp = @view height[1:len, 3:wid+2]
+    himjm = @view height[1:len, 1:wid]
+    hipjm = @view height[3:len+2, 1:wid]
+    
+    spreadingfactor = stats.γ * (1 - cospi(stats.θ)) * (stats.n-1)*(stats.m-1) / ((stats.n - stats.m)*stats.hstar)
+    
+    pressure = @strided -stats.γ .* ((2/3 .* (hip .+ him .+ hjp .+ hjm) .+ 
+                             1/6 .* (hipjp .+ himjp .+ himjm .+ hipjm) .- 
+                             20 .* h)) .- spreadingfactor .* (power_broad.((stats.hstar./(h .+ hᵪ)),stats.n) .- power_broad.((stats.hstar./(h .+ hᵪ)), stats.m)) 
+
+    return nothing
+end
 #= TODO: Test the GPU implementation 
 function capillarypressure!(pressure::CuArray, height, stats::pressurestats)
     # The default compute size    
@@ -307,10 +330,10 @@ function kernel_pressure!(p, height, linear, diagonal, stats::pressurestats, len
     him = view(height, Mod(0:len-1), 1:wid)
     hjp = view(height, 1:len, Mod(2:wid+1))
     hjm = view(height, 1:len, Mod(0:wid-1))
-    # hipjp = view(height, Mod(2:len+1), Mod(2:wid+1))
-    # himjp = view(height, Mod(0:len-1), Mod(2:wid+1))
-    # himjm = view(height, Mod(0:len-1), Mod(0:wid-1))
-    # hipjm = view(height, Mod(2:len+1), Mod(0:wid-1))
+    hipjp = view(height, Mod(2:len+1), Mod(2:wid+1))
+    himjp = view(height, Mod(0:len-1), Mod(2:wid+1))
+    himjm = view(height, Mod(0:len-1), Mod(0:wid-1))
+    hipjm = view(height, Mod(2:len+1), Mod(0:wid-1))
 
     linear .= (hip + him + hjp + hjm)
     # diagonal .= 1 ./ 6 .* (hipjp .+ himjp .+ himjm .+ hipjm)