Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GPU linear operators #163

Open
wants to merge 50 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
50 commits
Select commit Hold shift + click to select a range
5c2e71c
init design (no working implementation)
Abdelrahman912 Dec 18, 2024
8636a49
init setup
Abdelrahman912 Dec 20, 2024
34d8b79
make cuda work for the test kernel
Abdelrahman912 Dec 21, 2024
f7bd306
add gpuy subdof and init gpu update operator
Abdelrahman912 Jan 6, 2025
7eedd51
add subdof iterator
Abdelrahman912 Jan 6, 2025
16846e4
minor fix in sdh iterator
Abdelrahman912 Jan 6, 2025
5787e82
first working example (not refined)
Abdelrahman912 Jan 9, 2025
c58a14f
working example (not refined)
Abdelrahman912 Jan 9, 2025
324dbd0
Merge branch 'main' into gpu-operators
Abdelrahman912 Jan 9, 2025
139fc83
minor adjustment
Abdelrahman912 Jan 9, 2025
babc06c
minor refinement
Abdelrahman912 Jan 13, 2025
fe064cf
init add coeffs
Abdelrahman912 Jan 13, 2025
6bdae08
add to extension (doesn't dispatch to ext tho)
Abdelrahman912 Jan 14, 2025
c5c8b0b
minor edit
Abdelrahman912 Jan 14, 2025
4d3f7de
move to ext (working implementation)
Abdelrahman912 Jan 15, 2025
ac9c807
minor fix
Abdelrahman912 Jan 15, 2025
0d075a8
add coesfficients tests for gpu
Abdelrahman912 Jan 17, 2025
863ab05
minor refinements in coefficients
Abdelrahman912 Jan 17, 2025
0daf36c
minor fix
Abdelrahman912 Jan 17, 2025
7ba3139
restructure PR913
Abdelrahman912 Jan 21, 2025
077572f
change gpu -> device
Abdelrahman912 Jan 21, 2025
34c4536
some restructuring for memory allocation
Abdelrahman912 Jan 24, 2025
0e41a01
fix extension
Abdelrahman912 Jan 25, 2025
f512b88
optimize global mem
Abdelrahman912 Jan 25, 2025
cfc7cdd
minor fix for cuda op
Abdelrahman912 Jan 25, 2025
b59927d
pre kernel launch adapt for dh
Abdelrahman912 Feb 1, 2025
eb93c8d
Merge branch 'main' into gpu-operators
Abdelrahman912 Feb 1, 2025
d4d3bab
minor fix
Abdelrahman912 Feb 1, 2025
23f5413
init fix dh adapt error
Abdelrahman912 Feb 5, 2025
5e5216d
fix memory leak
Abdelrahman912 Feb 6, 2025
d50f7fa
pr review fix
Abdelrahman912 Feb 7, 2025
32c5763
init qp
Abdelrahman912 Feb 11, 2025
7bcdd88
minor fix
Abdelrahman912 Feb 11, 2025
8fe74d2
change operators and fuse mem alloc
Abdelrahman912 Feb 12, 2025
b6be5d4
subdof loop (i.e. get rid of unsafe cuda convert)
Abdelrahman912 Feb 13, 2025
110c7ca
Merge branch 'gpu-operators' into add-quadpoints-in-quadvalues
Abdelrahman912 Feb 13, 2025
ba95bae
use quad points instead of quad values
Abdelrahman912 Feb 13, 2025
97c6b00
minor fix
Abdelrahman912 Feb 14, 2025
69530cc
fix 1
Abdelrahman912 Feb 17, 2025
7339d1c
add complex test function and other stuff
Abdelrahman912 Feb 17, 2025
59fa5b4
fix CI (hopefully!)
Abdelrahman912 Feb 18, 2025
72a8f0e
minor fix
Abdelrahman912 Feb 18, 2025
569ef85
add n_threads and n_blocks as kwargs + change cuda default strategy n…
Abdelrahman912 Feb 18, 2025
a2a8bcd
fix SVector
Abdelrahman912 Feb 18, 2025
a0e28f4
add benchmarks + solve major bugs in memory allocation
Abdelrahman912 Feb 19, 2025
2215395
remove nested closures
Abdelrahman912 Feb 19, 2025
81a6529
add line + remove unnecessary code
Abdelrahman912 Feb 19, 2025
0c1d58f
Merge branch 'main' into gpu-operators
Abdelrahman912 Feb 19, 2025
d1bd648
add benchmarks
Abdelrahman912 Feb 20, 2025
291f3be
change analytical function + add PEALinearOperator benchmark
Abdelrahman912 Feb 20, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
minor fix
  • Loading branch information
Abdelrahman912 committed Jan 15, 2025
commit ac9c80701356e1d73b2e4e98dc2ed71df77b1904
3 changes: 2 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,8 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"


[weakdeps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

[extensions]
CuThunderboltExt = "CUDA"
Expand Down
18 changes: 6 additions & 12 deletions ext/cuda/cuda_adapt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,14 +22,7 @@ function Adapt.adapt_structure(to, cysc::CartesianCoordinateSystemCache)
return CartesianCoordinateSystemCache(cs, cv)
end

function Adapt.adapt_structure(to, cv::CellValues)
fv = Adapt.adapt(to, StaticInterpolationValues(cv.fun_values))
gm = Adapt.adapt(to, StaticInterpolationValues(cv.geo_mapping))
n_quadoints = cv.qr.weights |> length
weights = Adapt.adapt(to, ntuple(i -> cv.qr.weights[i], n_quadoints))
return StaticCellValues(fv, gm, weights)
end

# TODO: not used in the current codebase
function _convert_subdofhandler_to_gpu(cell_dofs, cell_dof_soffset, sdh::SubDofHandler)
GPUSubDofHandler(
cell_dofs,
Expand All @@ -41,10 +34,11 @@ function _convert_subdofhandler_to_gpu(cell_dofs, cell_dof_soffset, sdh::SubDofH
)
end

# TODO: here or in ferrite-addons?
function Adapt.adapt_structure(to, dh::DofHandler{sdim}) where sdim
grid = adapt_structure(to, dh.grid)
# field_names = Tuple(sym for sym in dh.field_names)
#IndexType = eltype(dh.cell_dofs)\
#IndexType = eltype(dh.cell_dofs)
#IndexVectorType = CuVector{IndexType}
cell_dofs = adapt(to, dh.cell_dofs .|> (i -> convert(Int32,i)) |> cu) # currently you cant create Dofhandler with Int32
cell_dofs_offset = adapt(to, dh.cell_dofs_offset .|> (i -> convert(Int32,i)) |> cu)
Expand All @@ -58,17 +52,17 @@ function Adapt.adapt_structure(to, dh::DofHandler{sdim}) where sdim
cell_dofs,
cell_dofs_offset,
cell_to_sdh,
dh.ndofs,
convert(Int32,dh.ndofs),
)
#return GPUDofHandler(dh, gpudata)
return GPUDofHandler(gpudata)
end



# TODO: here or in ferrite-addons?
function Adapt.adapt_structure(to, grid::Grid{sdim, cell_type, T}) where {sdim, cell_type, T}
node_type = typeof(first(grid.nodes))
cells = Adapt.adapt_structure(to, grid.cells |> cu)
cells = Adapt.adapt_structure(to, grid.cells .|> (x -> Int32.(x.nodes)) .|> eltype(grid.cells) |> cu)
nodes = Adapt.adapt_structure(to, grid.nodes |> cu)
#TODO subdomain info
return GPUGrid{sdim, cell_type, T, typeof(cells), typeof(nodes)}(cells, nodes)
Expand Down
12 changes: 3 additions & 9 deletions ext/cuda/cuda_operator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,10 @@ end

function _init_linop_cuda(linop::LinearOperator)
@unpack dh = linop
n_cells = dh |> get_grid |> getncells |> Int32
n_cells = dh |> get_grid |> getncells |> (x -> convert(Int32, x))
threads = convert(Int32, min(n_cells, 256))
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we need this to be In32?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe we don't need here Int32, I will fix it right away

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

actually I need them to be Int32, for allocating shared memory

function try_allocate_shared_mem(::Type{RHSObject{Tv}}, block_dim::Ti, n_basefuncs::Ti) where {Ti <: Integer, Tv <: Real}
    shared_mem = convert(Ti, sizeof(Tv) * (block_dim) * n_basefuncs)
    _can_use_dynshmem(shared_mem) || return nothing
    fe = DynamicSharedMemFunction{2, Tv, Ti}((block_dim, n_basefuncs), convert(Ti, 0))
    return SharedRHSMemAlloc( fe, shared_mem)
end

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this should rather be propagated by the assembly strategy type and not by the type of the input arguments.

blocks = _calculate_nblocks(threads, n_cells)
n_basefuncs = ndofs_per_cell(dh) |> Int32
n_basefuncs = convert(Int32,ndofs_per_cell(dh))
eles_caches = _setup_caches(linop)
mem_alloc = try_allocate_shared_mem(RHSObject{Float32}, threads, n_basefuncs)
mem_alloc isa Nothing || return CudaOperatorKernel(linop, threads, blocks, mem_alloc,eles_caches)
Expand Down Expand Up @@ -91,12 +91,11 @@ function Thunderbolt.update_operator!(op_ker::CudaOperatorKernel, time)
end



function _update_linear_operator_kernel!(b, dh_, eles_caches,mem_alloc, time)
dh = dh_.gpudata
for sdh_idx in 1:length(dh.subdofhandlers)
element_cache = eles_caches[sdh_idx]
for cell in CellIterator(dh,convert(Int32, sdh_idx) ,mem_alloc)
for cell in CellIterator(dh, convert(Int32,sdh_idx) ,mem_alloc)
bβ‚‘ = cellfe(cell)
assemble_element!(bβ‚‘, cell, element_cache, time)
dofs = celldofs(cell)
Expand All @@ -107,8 +106,3 @@ function _update_linear_operator_kernel!(b, dh_, eles_caches,mem_alloc, time)
end
return nothing
end


function Thunderbolt.test_ext(x::Float64)
println("Hello from the CUDA backend")
end
3 changes: 1 addition & 2 deletions src/Thunderbolt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,5 @@ export
PressureFieldBC,
BendingSpringBC,
RobinBC,
ConstantPressureBC,
test_ext
ConstantPressureBC
end
10 changes: 5 additions & 5 deletions src/ferrite-addons/PR913.jl
Original file line number Diff line number Diff line change
Expand Up @@ -304,11 +304,11 @@ end

Ferrite.CellIterator(dh::GPUDofHandlerData, buffer_alloc::AbstractGlobalMemAlloc) = _cell_iterator(dh, -1, dh |> get_grid |> getncells |> Int32, buffer_alloc) ## iterate over all cells

function Ferrite.CellIterator(dh::GPUDofHandlerData,sdh_idx::Integer, buffer_alloc::AbstractGlobalMemAlloc)
function Ferrite.CellIterator(dh::GPUDofHandlerData,sdh_idx::Ti, buffer_alloc::AbstractGlobalMemAlloc) where {Ti <: Integer}
## iterate over all cells in the subdomain
# check if the subdomain index is valid
sdh_idx βˆ‰ 1:length(dh.subdofhandlers) && return CudaOutOfBoundCellIterator()
n_cells = dh.subdofhandlers[sdh_idx].cellset |> length |> Int32
n_cells = dh.subdofhandlers[sdh_idx].cellset |> length |> (x -> convert(Ti, x))
return _cell_iterator(dh, sdh_idx,n_cells, buffer_alloc)
end

Expand All @@ -319,7 +319,7 @@ function Ferrite.CellIterator(dh::GPUDofHandlerData,sdh_idx::Integer, buffer_all
## iterate over all cells in the subdomain
# check if the subdomain index is valid
sdh_idx βˆ‰ 1:length(dh.subdofhandlers) && return CudaOutOfBoundCellIterator()
n_cells = dh.subdofhandlers[sdh_idx].cellset |> length |> Int32
n_cells = dh.subdofhandlers[sdh_idx].cellset |> length |> (x -> convert(typeof(dh.ndofs), x))
return _cell_iterator(dh, sdh_idx,n_cells, buffer_alloc)
end

Expand Down Expand Up @@ -353,7 +353,7 @@ struct GPUCellCache{Ti <: Integer, DOFS <: AbstractVector{Ti}, NN, NODES <: SVec
end


function _makecache(iterator::AbstractCUDACellIterator, e::Integer)
function _makecache(iterator::AbstractCUDACellIterator, e::Ti) where {Ti <: Integer}
dh = iterator.dh
grid = iterator.grid
sdh_idx = iterator.sdh_idx
Expand All @@ -362,7 +362,7 @@ function _makecache(iterator::AbstractCUDACellIterator, e::Integer)
cell = Ferrite.getcells(grid, e)

# Extract the node IDs of the cell.
nodes = SVector(convert.(Int32, Ferrite.get_node_ids(cell))...)
nodes = SVector(convert.(Ti, Ferrite.get_node_ids(cell))...)

# Extract the degrees of freedom for the cell.
dofs = Ferrite.celldofs(dh, e)
Expand Down
4 changes: 0 additions & 4 deletions src/gpu/gpu_operator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,3 @@ function update_operator!(::AbstractOperatorKernel, time)
error("Not implemented")
end


function test_ext(x::AbstractFloat)
error("Not implemented")
end
11 changes: 5 additions & 6 deletions test/gpu/operators-test.jl
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
## TODO: Put test operators here or with cpu operators?
using Thunderbolt
using CUDA
using Test
using StaticArrays



grid = generate_grid(Quadrilateral, (2,2))
left = Tensor{1, 2, Float32}((-1.0, -1.0)) # define the left bottom corner of the grid.
right = Tensor{1, 2, Float32}((1.0, 1.0)) # define the right top corner of the grid.
grid = generate_grid(Quadrilateral, (2,2),left,right)
dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefQuadrilateral,1}())
close!(dh)
Expand All @@ -19,8 +19,8 @@ propertynames(dh)


protocol = AnalyticalTransmembraneStimulationProtocol(
AnalyticalCoefficient((x,t) -> 1.0, CoordinateSystemCoefficient(cs)),
[SVector((0.0, 1.0))]
AnalyticalCoefficient((x,t) -> 1.f0, CoordinateSystemCoefficient(cs)),
[SVector((0.f0, 1.f0))]
)


Expand All @@ -43,7 +43,6 @@ cuda_op = Thunderbolt.init_linear_operator(CUDABackend,protocol, qrc, dh);
Thunderbolt.update_operator!(cuda_op,0.0)



@test Vector(cuda_op.op.b) β‰ˆ linop.b