Add type GPUSystem to determine if a system is stored on the CPU or GPU

examples/n_body/n_body_system.jl

@@ -1,7 +1,7 @@
 using TrixiParticles
 using LinearAlgebra
 
-struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS}
+struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothing}
     initial_condition :: InitialCondition{ELTYPE}
     mass              :: Array{ELTYPE, 1} # [particle]
     G                 :: ELTYPE

src/general/general.jl

@@ -1,12 +1,17 @@
-abstract type System{NDIMS} end
+# Abstract supertype for all system types. We additionally store the type of the system's
+# initial condition, which is `Nothing` when using KernelAbstractions.jl.
+abstract type System{NDIMS, IC} end
 
-abstract type FluidSystem{NDIMS} <: System{NDIMS} end
+# When using KernelAbstractions.jl, the initial condition has been replaced by `nothing`
+GPUSystem = System{NDIMS, Nothing} where {NDIMS}
+
+abstract type FluidSystem{NDIMS, IC} <: System{NDIMS, IC} end
 timer_name(::FluidSystem) = "fluid"
 
-abstract type SolidSystem{NDIMS} <: System{NDIMS} end
+abstract type SolidSystem{NDIMS, IC} <: System{NDIMS, IC} end
 timer_name(::SolidSystem) = "solid"
 
-abstract type BoundarySystem{NDIMS} <: System{NDIMS} end
+abstract type BoundarySystem{NDIMS, IC} <: System{NDIMS, IC} end
 timer_name(::BoundarySystem) = "boundary"
 
 @inline function set_zero!(du)

src/general/semidiscretization.jl

@@ -60,6 +60,8 @@ struct Semidiscretization{S, RU, RV, NS}
     end
 end
 
+GPUSemidiscretization = Semidiscretization{<:NTuple{<:Any, GPUSystem}}
+
 function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch,
                             periodic_box_min_corner=nothing,
                             periodic_box_max_corner=nothing, threaded_nhs_update=true)

src/schemes/boundary/system.jl

@@ -13,8 +13,7 @@ The interaction between fluid and boundary particles is specified by the boundar
 - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface.
    Note: currently it is assumed that all fluids have the same adhesion coefficient.
 """
-struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, CA} <:
-       BoundarySystem{NDIMS}
+struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem{NDIMS, IC}
     initial_condition    :: IC
     coordinates          :: CO # Array{ELTYPE, 2}
     boundary_model       :: BM
@@ -68,20 +67,22 @@ The interaction between fluid and boundary particles is specified by the boundar
     This is an experimental feature and may change in a future releases.
 
 """
-struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D} <: BoundarySystem{NDIMS}
-    coordinates      :: ARRAY2D # [dimension, particle]
-    radius           :: ARRAY1D # [particle]
-    normal_stiffness :: ELTYPE
+struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, IC,
+                         ARRAY1D, ARRAY2D} <: BoundarySystem{NDIMS, IC}
+    initial_condition :: IC
+    coordinates       :: ARRAY2D # [dimension, particle]
+    radius            :: ARRAY1D # [particle]
+    normal_stiffness  :: ELTYPE
 
     function BoundaryDEMSystem(initial_condition, normal_stiffness)
         coordinates = initial_condition.coordinates
         radius = 0.5 * initial_condition.particle_spacing *
                  ones(length(initial_condition.mass))
         NDIMS = size(coordinates, 1)
 
-        return new{NDIMS, eltype(coordinates), typeof(radius), typeof(coordinates)}(coordinates,
-                                                                                    radius,
-                                                                                    normal_stiffness)
+        return new{NDIMS, eltype(coordinates), typeof(initial_condition),
+                   typeof(radius), typeof(coordinates)}(initial_condition, coordinates,
+                                                        radius, normal_stiffness)
     end
 end
 

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -40,7 +40,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                     gravity-like source terms.
 """
 struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V,
-                                   PF, ST, C} <: FluidSystem{NDIMS}
+                                   PF, ST, C} <: FluidSystem{NDIMS, IC}
     initial_condition                 :: IC
     mass                              :: M # Vector{ELTYPE}: [particle]
     density_calculator                :: DC

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -42,7 +42,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method.
 
 """
 struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K,
-                                   V, DD, COR, PF, ST, SRFT, C} <: FluidSystem{NDIMS}
+                                   V, DD, COR, PF, ST, SRFT, C} <: FluidSystem{NDIMS, IC}
     initial_condition                 :: IC
     mass                              :: MA     # Array{ELTYPE, 1}
     pressure                          :: P      # Array{ELTYPE, 1}

src/schemes/solid/discrete_element_method/system.jl

@@ -26,8 +26,8 @@ specified material properties and contact mechanics.
  !!! warning "Experimental Implementation"
     This is an experimental feature and may change in a future releases.
 """
-struct DEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ST} <: SolidSystem{NDIMS}
-    initial_condition   :: InitialCondition{ELTYPE}
+struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST} <: SolidSystem{NDIMS, IC}
+    initial_condition   :: IC
     mass                :: ARRAY1D               # [particle]
     radius              :: ARRAY1D               # [particle]
     elastic_modulus     :: ELTYPE
@@ -53,7 +53,7 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ST} <: SolidSystem{NDIMS}
             throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
         end
 
-        return new{NDIMS, ELTYPE, typeof(mass),
+        return new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass),
                    typeof(source_terms)}(initial_condition, mass, radius, elastic_modulus,
                                          poissons_ratio, normal_stiffness,
                                          damping_coefficient, acceleration_, source_terms)

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -49,7 +49,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
     where `beam` and `fixed_particles` are of type `InitialCondition`.
 """
 struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D,
-                                K, PF, ST} <: SolidSystem{NDIMS}
+                                K, PF, ST} <: SolidSystem{NDIMS, IC}
     initial_condition   :: IC
     initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle]
     current_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle]

test/general/semidiscretization.jl

@@ -1,8 +1,8 @@
 # Use `@trixi_testset` to isolate the mock functions in a separate namespace
 @trixi_testset "Semidiscretization" begin
     # Mock systems
-    struct System1 <: TrixiParticles.System{3} end
-    struct System2 <: TrixiParticles.System{3} end
+    struct System1 <: TrixiParticles.System{3, Nothing} end
+    struct System2 <: TrixiParticles.System{3, Nothing} end
 
     system1 = System1()
     system2 = System2()
@@ -39,7 +39,7 @@
             struct BoundaryModelMock end
 
             # Mock fluid system
-            struct FluidSystemMock <: TrixiParticles.FluidSystem{2} end
+            struct FluidSystemMock <: TrixiParticles.FluidSystem{2, Nothing} end
 
             kernel = Val(:smoothing_kernel)
             Base.ndims(::Val{:smoothing_kernel}) = 2