Boundary conditions kernel and dependencies

.gitignore

@@ -4,3 +4,4 @@
 *.mp4
 Manifest.toml
 /docs/build/
+.vscode/

CUDAEnv/bcs.jl

@@ -0,0 +1,73 @@
+using KernelAbstractions, Adapt, OffsetArrays, BenchmarkTools
+
+if Base.find_package("CUDA") !== nothing
+    using CUDA
+    using CUDA.CUDAKernels
+    const backend = CUDABackend()
+    CUDA.allowscalar(false)
+else
+    const backend = CPU()
+end
+
+# helper functions
+ArrayT = (backend == CPU()) ? Array : CuArray
+splitn(n) = Base.front(n), n[end]
+size_u(u) = splitn(size(u))
+@inline δ(i, N::Int) = CartesianIndex(ntuple(j -> j == i ? 1 : 0, N))
+@inline δ(i, I::CartesianIndex{N}) where {N} = δ(i, N)
+@inline ∂(a, I::CartesianIndex{d}, f::AbstractArray{T,d}) where {T, d} = @inbounds f[I] - f[I - δ(a, I)]
+@inline ∂(a, I::CartesianIndex{m}, u::AbstractArray{T,n}) where {T, n, m} = @inbounds u[I + δ(a, I), a] - u[I, a]
+origin0(D=0) = OffsetArrays.Origin(D > 0 ? (zeros(Int, D)..., 1) : 0)
+function slice(dims::NTuple{N}, i, j, low = 1, trim = 0) where N
+    CartesianIndices(ntuple(k-> k == j ? (i:i) : (low:dims[k] - trim), N))
+end
+
+# flow and kernels
+struct Flow{B <: Backend, V, P, UBC, BC}
+    backend :: B
+    u :: V
+    σ :: P
+    U :: UBC
+    bcs :: BC
+end
+
+function Flow(N::NTuple{D}, U; backend = CPU(), ftype = Float32) where D
+    u = rand(ftype, ((N .+ 2)..., D)) |> origin0(D)
+    σ = Array{ftype}(undef, N .+ 2) |> origin0()
+    bcs, Ng = [], size(σ)
+    for d ∈ 1:D
+        slice_ghost_start = slice(Ng, 0, d, 1, 2)
+        slice_donor_start = slice_ghost_start .+ δ(d, D)
+        slice_ghost_end = slice(Ng, Ng[d] - 1, d, 1, 2)
+        slice_donor_end = slice_ghost_end .- δ(d, D)
+        push!(bcs, zip(slice_ghost_start, slice_donor_start, ntuple(x -> d, length(slice_ghost_start)))...,
+            zip(slice_ghost_end, slice_donor_end, ntuple(x -> d, length(slice_ghost_end)))...)
+    end
+    return Flow(backend, adapt(ArrayT, u), adapt(ArrayT, σ), U, adapt(ArrayT, Tuple.(bcs)))
+end
+
+# wrapper for boundary conditions kernel
+function bcs!(flow, f = 1.0)
+    _bcs!(backend, 64)(flow.u, flow.U, f, flow.bcs, ndrange=size(flow.bcs))
+end
+
+# bounary conditions kernel
+@kernel function _bcs!(u, @Const(U), @Const(f), @Const(bcs))
+    i = @index(Global, Linear)
+    ghostI, donorI, di = bcs[i][1], bcs[i][2], bcs[i][3]
+    _, D = size_u(u)
+    for d ∈ 1:D
+        if d == di
+            u[ghostI, d] = f * U[d]
+        else
+            u[ghostI, d] = u[donorI, d]
+        end
+    end
+end
+
+# main
+const FT = Float32
+N = (3, 4)
+
+flow = Flow(N, (1.0, 0.0, 0.0); backend = backend, ftype = FT)
+@btime bcs!($flow)

Project.toml

@@ -4,9 +4,13 @@ authors = ["Gabriel Weymouth <[email protected]>"]
 version = "0.2.4"
 
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"