Simple test for 3D Linear Advection Equation

examples/advection_basic_3d.jl

@@ -0,0 +1,39 @@
+# 3D Linear Advection Equation
+using Trixi, TrixiGPU
+using OrdinaryDiffEq
+
+advection_velocity = (0.2f0, -0.7f0, 0.5f0)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0f0, -1.0f0, -1.0f0)
+coordinates_max = (1.0f0, 1.0f0, 1.0f0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+tspan = (0.0f0, 1.0f0)
+
+# FIXME: Remember to export
+ode = TrixiGPU.semidiscretize_gpu(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, save_everystep = false, callback = callbacks);
+summary_callback()

src/solvers/dg_1d.jl

@@ -44,7 +44,8 @@ function weak_form_kernel!(du, derivative_dhat, flux_arr)
 end
 
 # Kernel for prolonging two interfaces
-function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids)
+function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids,
+                                    equations::AbstractEquations{1})
     j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
 
@@ -82,7 +83,8 @@ function surface_flux_kernel!(surface_flux_arr, interfaces_u,
 end
 
 # Kernel for setting interface fluxes
-function interface_flux_kernel!(surface_flux_values, surface_flux_arr, neighbor_ids)
+function interface_flux_kernel!(surface_flux_values, surface_flux_arr, neighbor_ids,
+                                equations::AbstractEquations{1})
     i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
 
@@ -153,16 +155,17 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms,
 end
 
 # Pack kernels for prolonging solution to interfaces
-function cuda_prolong2interfaces!(u, mesh::TreeMesh{1}, cache)
+function cuda_prolong2interfaces!(u, mesh::TreeMesh{1}, equations, cache)
     neighbor_ids = CuArray{Int64}(cache.interfaces.neighbor_ids)
     interfaces_u = CuArray{Float32}(cache.interfaces.u)
 
     size_arr = CuArray{Float32}(undef, size(interfaces_u, 2), size(interfaces_u, 3))
 
     prolong_interfaces_kernel = @cuda launch=false prolong_interfaces_kernel!(interfaces_u,
                                                                               u,
-                                                                              neighbor_ids)
-    prolong_interfaces_kernel(interfaces_u, u, neighbor_ids;
+                                                                              neighbor_ids,
+                                                                              equations)
+    prolong_interfaces_kernel(interfaces_u, u, neighbor_ids, equations;
                               configurator_2d(prolong_interfaces_kernel, size_arr)...,)
 
     cache.interfaces.u = interfaces_u  # copy back to host automatically
@@ -192,8 +195,9 @@ function cuda_interface_flux!(mesh::TreeMesh{1}, nonconservative_terms::False, e
 
     interface_flux_kernel = @cuda launch=false interface_flux_kernel!(surface_flux_values,
                                                                       surface_flux_arr,
-                                                                      neighbor_ids)
-    interface_flux_kernel(surface_flux_values, surface_flux_arr, neighbor_ids;
+                                                                      neighbor_ids,
+                                                                      equations)
+    interface_flux_kernel(surface_flux_values, surface_flux_arr, neighbor_ids, equations;
                           configurator_2d(interface_flux_kernel, size_arr)...,)
 
     cache.elements.surface_flux_values = surface_flux_values # copy back to host automatically
@@ -216,7 +220,7 @@ end
 
 # Pack kernels for calculating surface integrals
 function cuda_surface_integral!(du, mesh::TreeMesh{1}, equations, dg::DGSEM, cache)
-    # FIXME: `surface_integral`
+    # FIXME: Check `surface_integral`
     factor_arr = CuArray{Float32}([
                                       dg.basis.boundary_interpolation[1, 1],
                                       dg.basis.boundary_interpolation[size(du, 2), 2]
@@ -262,7 +266,7 @@ function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{1}, equations, initial_condit
     cuda_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations,
                           dg.volume_integral, dg)
 
-    cuda_prolong2interfaces!(u, mesh, cache)
+    cuda_prolong2interfaces!(u, mesh, equations, cache)
 
     cuda_interface_flux!(mesh, have_nonconservative_terms(equations), equations, dg, cache)
 

src/solvers/dg_2d.jl

@@ -4,7 +4,7 @@
 # the @cuda macro with parameters from the kernel configurator. They are purely run on 
 # the device (i.e., GPU).
 
-# Kernel for calculating fluxes along normal direction
+# Kernel for calculating fluxes along normal directions
 function flux_kernel!(flux_arr1, flux_arr2, u, equations::AbstractEquations{2},
                       flux::Function)
     j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
@@ -52,7 +52,8 @@ function weak_form_kernel!(du, derivative_dhat, flux_arr1, flux_arr2)
 end
 
 # Kernel for prolonging two interfaces 
-function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids, orientations)
+function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids, orientations,
+                                    euqations::AbstractEquations{2})
     j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
 
@@ -104,7 +105,7 @@ end
 
 # Kernel for setting interface fluxes
 function interface_flux_kernel!(surface_flux_values, surface_flux_arr, neighbor_ids,
-                                orientations)
+                                orientations, equations::AbstractEquations{2})
     i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
     k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
@@ -198,7 +199,7 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{2}, nonconservative_terms,
 end
 
 # Pack kernels for prolonging solution to interfaces
-function cuda_prolong2interfaces!(u, mesh::TreeMesh{2}, cache)
+function cuda_prolong2interfaces!(u, mesh::TreeMesh{2}, equations, cache)
     neighbor_ids = CuArray{Int64}(cache.interfaces.neighbor_ids)
     orientations = CuArray{Int64}(cache.interfaces.orientations)
     interfaces_u = CuArray{Float32}(cache.interfaces.u)
@@ -210,8 +211,9 @@ function cuda_prolong2interfaces!(u, mesh::TreeMesh{2}, cache)
     prolong_interfaces_kernel = @cuda launch=false prolong_interfaces_kernel!(interfaces_u,
                                                                               u,
                                                                               neighbor_ids,
-                                                                              orientations)
-    prolong_interfaces_kernel(interfaces_u, u, neighbor_ids, orientations;
+                                                                              orientations,
+                                                                              euqations)
+    prolong_interfaces_kernel(interfaces_u, u, neighbor_ids, orientations, euqations;
                               configurator_2d(prolong_interfaces_kernel, size_arr)...,)
 
     cache.interfaces.u = interfaces_u  # copy back to host automatically
@@ -245,9 +247,10 @@ function cuda_interface_flux!(mesh::TreeMesh{2}, nonconservative_terms::False, e
     interface_flux_kernel = @cuda launch=false interface_flux_kernel!(surface_flux_values,
                                                                       surface_flux_arr,
                                                                       neighbor_ids,
-                                                                      orientations)
-    interface_flux_kernel(surface_flux_values, surface_flux_arr, neighbor_ids, orientations;
-                          configurator_3d(interface_flux_kernel, size_arr)...,)
+                                                                      orientations,
+                                                                      equations)
+    interface_flux_kernel(surface_flux_values, surface_flux_arr, neighbor_ids, orientations,
+                          equations; configurator_3d(interface_flux_kernel, size_arr)...,)
 
     cache.elements.surface_flux_values = surface_flux_values # copy back to host automatically
 
@@ -317,7 +320,7 @@ function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{2}, equations, initial_condit
     cuda_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations,
                           dg.volume_integral, dg)
 
-    cuda_prolong2interfaces!(u, mesh, cache)
+    cuda_prolong2interfaces!(u, mesh, equations, cache)
 
     cuda_interface_flux!(mesh, have_nonconservative_terms(equations), equations, dg, cache)