Format

src/TrixiGPU.jl

@@ -3,7 +3,8 @@ module TrixiGPU
 # Include other packages that are used in TrixiGPU.jl
 # using Reexport: @reexport
 
-using CUDA: @cuda, CuArray, HostKernel, launch_configuration, threadIdx
+using CUDA: @cuda, CuArray, HostKernel, launch_configuration,
+            threadIdx, blockIdx, blockDim
 using Trixi: AbstractEquations
 
 import Trixi: get_node_vars, get_node_coords, get_surface_node_vars
@@ -13,6 +14,7 @@ using StrideArrays: PtrArray
 # Include other source files
 include("function.jl")
 include("auxiliary/auxiliary.jl")
+include("solvers/solvers.jl")
 
 # Export the public APIs
 # export configurator_1d, configurator_2d, configurator_3d

src/solvers/dg_1d.jl

@@ -1,2 +1,60 @@
 # Solver functions for 1D DG methods
 
+function flux_kernel!(flux_arr, u, flux::Function, equations::AbstractEquations{1})
+    j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+
+    if (j <= size(u, 2) && k <= size(u, 3))
+        u_node = get_nodes_vars(u, equations, j, k)
+
+        flux_node = flux(u_node, 1, equations)
+
+        @inbounds begin
+            for ii in axes(u, 1)
+                flux_arr[ii, j, k] = flux_node[ii]
+            end
+        end
+    end
+
+    return nothing
+end
+
+function weak_form_kernel!(du, derivative_dhat, flux_arr, equations::AbstractEquations{1})
+    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+    k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
+
+    if (i <= size(du, 1) && j <= size(du, 2) && k <= size(du, 3))
+        @inbounds begin
+            for ii in axes(du, 2)
+                du[i, j, k] += derivative_dhat[j, ii] * flux_arr[i, ii, k]
+            end
+        end
+    end
+
+    return nothing
+end
+
+function volume_flux_kernel!(volume_flux_arr, u, volume_flux::Function,
+                             equations::AbstractEquations{1})
+    j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+
+    if (j <= size(u, 2)^2 && k <= size(u, 3))
+        j1 = div(j - 1, size(u, 2)) + 1
+        j2 = rem(j - 1, size(u, 2)) + 1
+
+        u_node = get_nodes_vars(u, equations, j1, k)
+        u_node1 = get_nodes_vars(u, equations, j2, k)
+
+        volume_flux_node = volume_flux(u_node, u_node1, 1, equations)
+
+        @inbounds begin
+            for ii in axes(u, 1)
+                volume_flux_arr[ii, j1, j2, k] = volume_flux_node[ii]
+            end
+        end
+    end
+
+    return nothing
+end

src/solvers/solvers.jl

@@ -1 +1,4 @@
 include("common.jl")
+include("dg_1d.jl")
+include("dg_2d.jl")
+include("dg_3d.jl")

test/test_solvers.jl

@@ -0,0 +1,22 @@
+using Trixi, TrixiGPU
+using CUDA: @cuda, CuArray
+using Test
+
+# @testset "Test solver functions" begin
+
+#     weak_form_kernel = @cuda launch=false weak_form_kernel!(du, derivative_dhat, flux_arr)
+# end
+
+advection_velocity = 1.0f0
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+du = CuArray{Float64}(undef, 10, 10, 10)
+derivative_dhat = CuArray{Float64}(undef, 10, 10)
+flux_arr = CuArray{Float64}(undef, 10, 10, 10)
+
+weak_form_kernel = @cuda launch=false TrixiGPU.weak_form_kernel!(du, derivative_dhat,
+                                                                 flux_arr, equations)
+weak_form_kernel(du,
+                 derivative_dhat,
+                 flux_arr, equations;
+                 TrixiGPU.configurator_3d(weak_form_kernel, du)...,)