Skip to content

Commit

Permalink
Add GPU tests and fix Compton implementation for the testss
Browse files Browse the repository at this point in the history
  • Loading branch information
AntonReinhard committed May 27, 2024
1 parent c768a57 commit fb8dee1
Show file tree
Hide file tree
Showing 12 changed files with 226 additions and 42 deletions.
13 changes: 5 additions & 8 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,15 +1,10 @@
name = "QEDprocesses"
uuid = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
authors = [
"Uwe Hernandez Acosta <u.hernandez@hzdr.de>",
"Simeon Ehrig",
"Klaus Steiniger",
"Tom Jungnickel",
"Anton Reinhard",
]
authors = ["Uwe Hernandez Acosta <u.hernandez@hzdr.de>", "Simeon Ehrig", "Klaus Steiniger", "Tom Jungnickel", "Anton Reinhard"]
version = "0.1.0"

[deps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Expand All @@ -21,10 +16,12 @@ StaticArrays = "1"
julia = "1.6"

[extras]
AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Random", "SafeTestsets", "Suppressor", "Test"]
test = ["Random", "SafeTestsets", "Suppressor", "Test", "AMDGPU", "CUDA"]
1 change: 1 addition & 0 deletions src/QEDprocesses.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ export AbstractModelDefinition, fundamental_interaction_type
export AbstractProcessDefinition, incoming_particles, outgoing_particles
export number_incoming_particles, number_outgoing_particles
export particles, number_particles
export in_phase_space_dimension, out_phase_space_dimension

# probabilities
export differential_probability, unsafe_differential_probability
Expand Down
12 changes: 6 additions & 6 deletions src/momentum_generation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,12 @@ function _generate_outgoing_momenta end

"""
_generate_momenta(
proc::AbstractProcessDefinition,
model::AbstractModelDefinition,
phase_space_def::AbstractPhasespaceDefinition,
in_phase_space::NTuple{N,T},
out_phase_space::NTuple{M,T},
) where {N,M,T<:Real}
proc::AbstractProcessDefinition,
model::AbstractModelDefinition,
phase_space_def::AbstractPhasespaceDefinition,
in_phase_space::NTuple{N,T},
out_phase_space::NTuple{M,T},
) where {N,M,T<:Real}
Return four-momenta for incoming and outgoing particles for given coordinate based phase-space points.
"""
Expand Down
10 changes: 5 additions & 5 deletions src/patch_QEDbase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
#
#############

# fix: https://github.com/QEDjl-project/QEDbase.jl/pull/61
# fix: https://github.com/QEDjl-project/QEDbase.jl/pull/65
Base.show(io::IO, ::Electron) = print(io, "electron")
Base.show(io::IO, ::Positron) = print(io, "positron")
Base.show(io::IO, ::Photon) = print(io, "photon")
Expand All @@ -24,7 +24,7 @@ Broadcast.broadcastable(part::AbstractParticleType) = Ref(part)
Broadcast.broadcastable(spin_or_pol::AbstractSpinOrPolarization) = Ref(spin_or_pol)

# fix: https://github.com/QEDjl-project/QEDbase.jl/pull/63
number_of_spin_pol(::AbstractDefinitePolarization) = 1
number_of_spin_pol(::AbstractDefiniteSpin) = 1
number_of_spin_pol(::AbstractIndefinitePolarization) = 2
number_of_spin_pol(::AbstractIndefiniteSpin) = 2
multiplicity(::AbstractDefinitePolarization) = 1
multiplicity(::AbstractDefiniteSpin) = 1
multiplicity(::AbstractIndefinitePolarization) = 2
multiplicity(::AbstractIndefiniteSpin) = 2
12 changes: 10 additions & 2 deletions src/phase_spaces/access.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,16 @@ momentum(part::ParticleStateful) = part.mom
Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`.
"""
momenta(psp::PhaseSpacePoint, ::Incoming) = momentum.(psp.in_particles)
momenta(psp::PhaseSpacePoint, ::Outgoing) = momentum.(psp.out_particles)
function momenta(
psp::PhaseSpacePoint{P,M,D,I,O,E}, ::Incoming
) where {N,P,M,D,I<:Tuple{Vararg{ParticleStateful,N}},O,E}
return NTuple{N,E}(momentum(p) for p in psp.in_particles)
end
function momenta(
psp::PhaseSpacePoint{P,M,D,I,O,E}, ::Outgoing
) where {N,P,M,D,I,O<:Tuple{Vararg{ParticleStateful,N}},E}
return NTuple{N,E}(momentum(p) for p in psp.out_particles)
end

"""
Base.getindex(psp::PhaseSpacePoint, dir::Incoming, n::Int)
Expand Down
33 changes: 14 additions & 19 deletions src/processes/one_photon_compton/perturbative/cross_section.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ end
We average over the initial spins and pols, and sum over final.
"""
function _averaging_norm(proc::Compton)
normalizations = number_of_spin_pol.(_in_spin_and_pol(proc))
normalizations = multiplicity.(_in_spin_and_pol(proc))
return inv(prod(normalizations))
end

Expand Down Expand Up @@ -74,8 +74,8 @@ end
in_photon_state = base_state(Photon(), Incoming(), in_photon_mom, proc.in_pol)

out_electron_state = base_state(Electron(), Outgoing(), out_electron_mom, proc.out_spin)

out_photon_state = base_state(Photon(), Outgoing(), out_photon_mom, proc.out_pol)

return _pert_compton_matrix_element(
in_electron_mom,
in_electron_state,
Expand Down Expand Up @@ -105,23 +105,18 @@ function _pert_compton_matrix_element(
QEDbase._as_svec(out_photon_state),
)

matrix_elements = Vector{ComplexF64}()
sizehint!(matrix_elements, length(base_states_comb))
for (in_el, in_ph, out_el, out_ph) in base_states_comb
push!(
matrix_elements,
_pert_compton_matrix_element_single(
in_electron_mom,
in_el,
in_photon_mom,
in_ph,
out_electron_mom,
out_el,
out_photon_mom,
out_ph,
),
)
end
matrix_elements = NTuple{length(base_states_comb),ComplexF64}(
_pert_compton_matrix_element_single(
in_electron_mom,
in_el,
in_photon_mom,
in_ph,
out_electron_mom,
out_el,
out_photon_mom,
out_ph,
) for (in_el, in_ph, out_el, out_ph) in base_states_comb
)

return matrix_elements
end
Expand Down
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@

function _total_probability(in_psp::InPhaseSpacePoint{<:Compton,<:PerturbativeQED})
omega = getE(momentum(in_psp[Incoming(), 2]))
omega = getE(momentum(in_psp, Incoming(), 2))

function func(x)
return unsafe_differential_probability(
Expand All @@ -9,7 +9,6 @@ function _total_probability(in_psp::InPhaseSpacePoint{<:Compton,<:PerturbativeQE
end

tot_prob, _ = quadgk(func, -1, 1; rtol=sqrt(eps(omega)))

tot_prob *= 2 * pi # phi integration is trivial
return tot_prob
end
163 changes: 163 additions & 0 deletions test/gpu/process_interface.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,163 @@
using QEDprocesses
using QEDbase
using AMDGPU
using CUDA
using Random
using SafeTestsets

include("../test_implementation/TestImplementation.jl")

GPU_VECTOR_TYPES = Vector{Type}()
if AMDGPU.functional()
println("Testing with AMDGPU.jl")
push!(GPU_VECTOR_TYPES, ROCVector)
end
if CUDA.functional()
println("Testing with CUDA.jl")
push!(GPU_VECTOR_TYPES, CuVector)
end

if isempty(GPU_VECTOR_TYPES)
println("No functional GPUs found!")
return nothing
end

PROC_DEF_TUPLES = [
(
Compton(),
PerturbativeQED(),
PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
),
(
Compton(SpinUp(), PolX(), SpinDown(), PolY()),
PerturbativeQED(),
PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
),
]

RNG = Random.MersenneTwister(573)

@testset "Testing with $VECTOR_TYPE" for VECTOR_TYPE in GPU_VECTOR_TYPES
@testset "$proc $model $ps_def" for (proc, model, ps_def) in PROC_DEF_TUPLES
N = 10_000

psps = [
PhaseSpacePoint(
proc,
model,
ps_def,
TestImplementation._rand_coordinates(RNG, proc, model, ps_def)...,
) for _ in 1:N
]
procs = [proc for _ in 1:N]

gpupsps = VECTOR_TYPE(psps)
gpuprocs = VECTOR_TYPE(procs)

@testset "PSP interface" begin
in_moms_gpu = Vector(momenta.(gpupsps, Incoming()))
out_moms_gpu = Vector(momenta.(gpupsps, Outgoing()))
in_moms = momenta.(psps, Incoming())
out_moms = momenta.(psps, Outgoing())

@test getindex.(in_moms_gpu, Ref(1)) == getindex.(in_moms, Ref(1))
@test getindex.(in_moms_gpu, Ref(2)) == getindex.(in_moms, Ref(2))
@test getindex.(out_moms_gpu, Ref(1)) == getindex.(out_moms, Ref(1))
@test getindex.(out_moms_gpu, Ref(2)) == getindex.(out_moms, Ref(2))
end

@testset "Private Process Functions" begin
@test all(
isapprox.(
Vector(QEDprocesses._averaging_norm.(gpuprocs)),
QEDprocesses._averaging_norm.(procs),
),
)
end

@testset "Public Process Functions" begin
@test Vector(incoming_particles.(gpuprocs)) == incoming_particles.(procs)
@test Vector(outgoing_particles.(gpuprocs)) == outgoing_particles.(procs)

@test Vector(particles.(gpuprocs, Incoming())) == particles.(procs, Incoming())
@test Vector(particles.(gpuprocs, Outgoing())) == particles.(procs, Outgoing())

@test Vector(number_incoming_particles.(gpuprocs)) ==
number_incoming_particles.(procs)
@test Vector(number_outgoing_particles.(gpuprocs)) ==
number_outgoing_particles.(procs)

@test Vector(number_particles.(gpuprocs, Incoming())) ==
number_particles.(procs, Incoming())
@test Vector(number_particles.(gpuprocs, Outgoing())) ==
number_particles.(procs, Outgoing())

@test Vector(in_phase_space_dimension.(gpuprocs, model)) ==
in_phase_space_dimension.(procs, model)
@test Vector(out_phase_space_dimension.(gpuprocs, model)) ==
out_phase_space_dimension.(procs, model)
end

@testset "Private PhaseSpacePoint Interface" begin
@test all(
isapprox.(
Vector(QEDprocesses._incident_flux.(gpupsps)),
QEDprocesses._incident_flux.(psps),
),
)

@test all(
TestImplementation.tuple_isapprox.(
Vector(QEDprocesses._matrix_element.(gpupsps)),
QEDprocesses._matrix_element.(psps),
),
)

@test Vector(QEDprocesses._is_in_phasespace.(gpupsps)) ==
QEDprocesses._is_in_phasespace.(psps)

@test all(
isapprox.(
Vector(QEDprocesses._phase_space_factor.(gpupsps)),
QEDprocesses._phase_space_factor.(psps),
),
)

#= TODO: this is currently broken
@test all(
isapprox.(
Vector(QEDprocesses._total_probability.(gpupsps)),
QEDprocesses._total_probability.(psps),
),
)=#
end

@testset "Public PhaseSpacePoint Interface" begin
@test all(
isapprox.(
Vector(differential_probability.(gpupsps)),
differential_probability.(psps),
),
)

@test all(
isapprox.(
Vector(QEDprocesses._is_in_phasespace.(gpupsps)),
QEDprocesses._is_in_phasespace.(psps),
),
)

@test all(
isapprox.(
Vector(differential_cross_section.(gpupsps)),
differential_cross_section.(psps),
),
)

#= TODO: this is currently broken
@test all(
isapprox.(Vector(total_cross_section.(gpupsps)), total_cross_section.(psps))
)=#
end
end
end
4 changes: 4 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,4 +28,8 @@ begin

# scattering processes
include("processes/run_process_test.jl")

@time @safetestset "GPU process interface" begin
include("gpu/process_interface.jl")
end
end
1 change: 1 addition & 0 deletions test/test_implementation/TestImplementation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ include("groundtruths.jl")
include("test_model.jl")
include("test_process.jl")
include("random_momenta.jl")
include("random_coordinates.jl")
include("utils.jl")

end
8 changes: 8 additions & 0 deletions test/test_implementation/random_coordinates.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
"""
Return a tuple of tuples of incoming and outgoing coordinates for a given process, model and ps_def that make up a physical phase space point.
"""
function _rand_coordinates(
rng::AbstractRNG, ::PROCESS, ::MODEL, ::PS_DEF
) where {PROCESS<:Compton,MODEL<:PerturbativeQED,PS_DEF<:PhasespaceDefinition}
return ((rand(rng, Float64),), (rand(rng, Float64), rand(rng, Float64)))
end
8 changes: 8 additions & 0 deletions test/test_implementation/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,3 +42,11 @@ end
function _furl_moms(moms::NTuple{N,Float64}) where {N}
return Tuple(_furl_moms(Vector{Float64}([moms...])))
end

tuple_isapprox(::Tuple{}, ::Tuple{}; atol=0.0, rtol=eps()) = true
function tuple_isapprox(
a::Tuple{<:Number,Vararg}, b::Tuple{<:Number,Vararg}; atol=0.0, rtol=eps()
)
return isapprox(a[1], b[1]; atol=atol, rtol=rtol) &&
tuple_isapprox(a[2:end], b[2:end]; atol=atol, rtol=rtol)
end

0 comments on commit fb8dee1

Please sign in to comment.