From 93c65becf31badd2f02ae703cc89aa50d21319db Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sun, 11 Oct 2026 20:41:47 +0200 Subject: [PATCH 1/5] Give ParticleList a species table, a charge and a mass MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add Species, with a name, a charge and a mass, and give ParticleList a species table and a per-particle species index, with the accessors charge, mass and species. Add the lazy chargeweights matrix of q·w for a charge deposit, and the species group of the HDF5 round trip. A list built without a species holds one species with charge and mass one, so q·w is w and existing results are unchanged. --- CHANGELOG.md | 10 +- src/ParticleMethods.jl | 8 ++ src/particle_list.jl | 133 ++++++++++++++++++++++-- src/species.jl | 61 +++++++++++ test/particle_list.jl | 222 +++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/species.jl | 39 ++++++++ 7 files changed, 465 insertions(+), 9 deletions(-) create mode 100644 src/species.jl create mode 100644 test/species.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 8d40c73..e1b99f2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,14 @@ first entry is written. in-place update over `pl.x` and `pl.v` under `--track-allocation`. It measures rather than asserts, so `test/` was the wrong home: a profiling harness that runs as part of the suite costs a minute and checks nothing. +- **`ParticleList` carries a species table and a per-particle species index.** A new `Species` + type holds a name, a charge and a mass. `charge(pl, i)`, `mass(pl, i)` and `species(pl, i)` + read them back for each particle, and `chargeweights(pl)` returns the charge-weighted weights + `q_i · w_i` as a lazy `1 × np` matrix, the weight a charge deposit reads. A list built without + a species holds one species with charge and mass one, so `q · w` is `w` and every existing + result is unchanged. `h5save` writes the species table and the index as a group beside + `particles`, and a file without that group loads with the default species. The species index + costs 8 bytes per particle at construction. ### Changed @@ -48,7 +56,7 @@ first entry is written. per particle — a redundant representation of the underlying matrix. It now uses `eachslice` over each variable's view, yielding a zero-allocation `AbstractVector` with the same element type and values. Measured at `np = 10⁶` particles, `nd = 6` - components, three variables: allocation reduced from 328.0 MB to 184.0 MB (43.9% + components, three variables: allocation reduced from 328.0 MB to 192.0 MB (41.5% reduction). - **`StaticArrays` removed from `[deps]` in Project.toml.** The `using StaticArrays: MVector` import in `src/ParticleMethods.jl` was the only use in `src/` and was removed by the `Particle(DT, len)` diff --git a/src/ParticleMethods.jl b/src/ParticleMethods.jl index 67b86ca..c810aec 100644 --- a/src/ParticleMethods.jl +++ b/src/ParticleMethods.jl @@ -11,6 +11,13 @@ using HDF5: H5DataStore include("hdf5_utils.jl") +include("species.jl") + +export Species +export charge +export mass +export species + include("particle.jl") export Particle @@ -19,5 +26,6 @@ include("particle_list.jl") export ParticleList export eachparticle +export chargeweights end # module diff --git a/src/particle_list.jl b/src/particle_list.jl index 50ef860..f713226 100644 --- a/src/particle_list.jl +++ b/src/particle_list.jl @@ -8,7 +8,8 @@ _sort_names(::NamedTuple{N, T}) where {N, T} = Tuple(sort([N...])) _sort_ntuple(nt::NamedTuple) = NamedTuple{_sort_names(nt)}(nt) """ - ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple()) + ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(), + species = nothing, speciesindex = nothing) ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple()) ParticleList(x, v, w; parameters = NamedTuple()) ParticleList(h5::H5DataStore, path = "/") @@ -20,7 +21,14 @@ particle, and `pl[i]` is that particle as a [`Particle`](@ref). `variables` maps names to row indices or row ranges of `list`, and `parameters` maps names to values. Each name is a property of the list: `pl.x` is a view into the rows of a variable `x` for all particles, and the stored value for a parameter. The names must not repeat each other or the -field names `list`, `views`, `params`, `particles`, `variables` and `indices`. +field names `list`, `views`, `params`, `particles`, `variables`, `indices`, `species` and +`speciesindex`. + +`species` is the table of particle species, each a [`Species`](@ref) with a charge and a mass, and +`speciesindex` gives the species of each particle as an index into that table. A list built +without `species` holds one species with charge and mass one, so `chargeweights(pl)` is `pl.w`. A +list with more than one species needs an explicit `speciesindex` of length `np`, with every index +in `1:length(species)`. The second form makes a zero list of element type `DT` with `np` particles of `nd` components. The third form stacks positions `x`, velocities `v` and weights `w`, and defines the variables @@ -35,24 +43,51 @@ struct ParticleList{T, ST <: AbstractMatrix{T}, VT <: NamedTuple, PT <: NamedTup particles::PART variables::VART indices::IND + species::Vector{Species{T}} + speciesindex::Vector{Int} - function ParticleList(list, views, params, particles, variables, indices) + function ParticleList(list, views, params, particles, variables, indices, species, + speciesindex) @assert isdisjoint(fieldnames(ParticleList), keys(views)) @assert isdisjoint(fieldnames(ParticleList), keys(params)) @assert isdisjoint(keys(views), keys(params)) new{eltype(list), typeof(list), typeof(views), typeof(params), typeof(particles), typeof(variables), typeof(indices)}( - list, views, params, particles, variables, indices) + list, views, params, particles, variables, indices, species, speciesindex) end end -function ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple()) +function ParticleList( + list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(), + species = nothing, speciesindex = nothing) svariables = _sort_ntuple(variables) views = map(idx_range -> view(list, idx_range, :), svariables) vars = map(v -> eachslice(v; dims = ndims(v)), views) particles = [Particle(p; variables = svariables, parameters = parameters) for p in eachcol(list)] - ParticleList(list, views, parameters, particles, vars, svariables) + + T = eltype(list) + np = size(list, 2) + stable = species === nothing ? [Species(:default, one(T), one(T))] : + Species{T}.(species) + sind = if speciesindex === nothing + if length(stable) > 1 + throw(ArgumentError( + "a ParticleList with $(length(stable)) species needs an explicit speciesindex")) + end + ones(Int, np) + else + convert(Vector{Int}, speciesindex) + end + + length(sind) == np || throw(ArgumentError( + "speciesindex has length $(length(sind)), expected $np")) + for i in sind + 1 ≤ i ≤ length(stable) || throw(ArgumentError( + "speciesindex contains $i, outside 1:$(length(stable))")) + end + + ParticleList(list, views, parameters, particles, vars, svariables, stable, sind) end function ParticleList(DT::DataType, np::Int, nd::Int; kwargs...) @@ -84,7 +119,9 @@ function Base.:(==)(pl1::ParticleList{T1, ST1}, pl2::ParticleList{ && pl1.params == pl2.params && pl1.particles == pl2.particles && pl1.variables == pl2.variables - && pl1.indices == pl2.indices) + && pl1.indices == pl2.indices + && pl1.species == pl2.species + && pl1.speciesindex == pl2.speciesindex) end @inline function Base.hasproperty(::ParticleList{T, ST, VT, PT}, s::Symbol) where { @@ -135,6 +172,68 @@ eachparticle(pl::ParticleList) = pl.particles Base.iterate(pl::ParticleList, state...) = iterate(pl.particles, state...) +""" + species(pl::ParticleList, i::Integer) + +Return the [`Species`](@ref) of particle `i` in `pl`, that is `pl.species[pl.speciesindex[i]]`. +""" +@inline species(pl::ParticleList, i::Integer) = pl.species[pl.speciesindex[i]] + +""" + charge(pl::ParticleList, i::Integer) + +Return the charge of particle `i` in `pl`, the charge of its [`Species`](@ref). +""" +@inline charge(pl::ParticleList, i::Integer) = species(pl, i).charge + +""" + mass(pl::ParticleList, i::Integer) + +Return the mass of particle `i` in `pl`, the mass of its [`Species`](@ref). +""" +@inline mass(pl::ParticleList, i::Integer) = species(pl, i).mass + +""" + ChargeWeights{T} <: AbstractMatrix{T} + +The charge-weighted weights `q_i · w_i` of the particles in a [`ParticleList`](@ref), as a lazy +`1 × np` matrix whose `[1, i]` element is the charge of particle `i` times its weight. Returned by +[`chargeweights`](@ref). +""" +struct ChargeWeights{T, PL <: ParticleList} <: AbstractMatrix{T} + parent::PL + row::Int +end + +Base.size(cw::ChargeWeights) = (1, length(cw.parent)) + +@inline function Base.getindex(cw::ChargeWeights, i::Int, j::Int) + @boundscheck checkbounds(cw, i, j) + return charge(cw.parent, j) * cw.parent.list[cw.row, j] +end + +""" + chargeweights(pl::ParticleList) + +Return the charge-weighted weights of the particles in `pl` as a lazy `1 × np` matrix, whose +`[1, i]` element is `charge(pl, i) * pl.w[1, i]`. This is the weight a charge deposit reads, so a +deposit that weights each particle by its charge changes one argument and no signature. + +It throws an `ArgumentError` when `pl` has no variable `w`, or when `w` spans more than one row. +""" +function chargeweights(pl::ParticleList{T}) where {T} + haskey(pl.indices, :w) || throw(ArgumentError("the ParticleList has no variable w")) + return ChargeWeights{T, typeof(pl)}(pl, _wrow(pl.indices.w)) +end + +_wrow(w::Integer) = Int(w) + +function _wrow(w) + length(w) == 1 || + throw(ArgumentError("the variable w must span exactly one row, got $(length(w))")) + return Int(first(w)) +end + function ParticleList(h5::H5DataStore, path::AbstractString = "/") group = h5[path] @@ -151,7 +250,19 @@ function ParticleList(h5::H5DataStore, path::AbstractString = "/") params = NamedTuple{paramkeys}(paramvals) - ParticleList(particles; variables = variables, parameters = params) + species, speciesindex = if haskey(group, "species") + sgroup = group["species"] + names = Symbol.(read(sgroup["name"])) + charges = read(sgroup["charge"]) + masses = read(sgroup["mass"]) + ([Species(n, q, m) for (n, q, m) in zip(names, charges, masses)], + convert(Vector{Int}, read(sgroup["index"]))) + else + (nothing, nothing) + end + + ParticleList(particles; variables = variables, parameters = params, + species = species, speciesindex = speciesindex) end function ParticleList(fpath::AbstractString, path::AbstractString = "/") @@ -175,6 +286,12 @@ function h5save(h5::H5DataStore, p::ParticleList; path::AbstractString = "/") params[string(key)] = p.params[key] end + sgroup = _create_group(group, "species") + sgroup["name"] = String[String(s.name) for s in p.species] + sgroup["charge"] = [s.charge for s in p.species] + sgroup["mass"] = [s.mass for s in p.species] + sgroup["index"] = p.speciesindex + group["particles"] = p.list end diff --git a/src/species.jl b/src/species.jl new file mode 100644 index 0000000..0b89399 --- /dev/null +++ b/src/species.jl @@ -0,0 +1,61 @@ +""" + Species(name::Symbol, charge, mass) + Species{T}(name::Symbol, charge, mass) + +A particle species: a name, a charge `q` and a mass `m`, shared by every particle of that species. +`Species(name, charge, mass)` promotes the charge and the mass to their common type, and +`Species{T}(s::Species)` converts an existing species to the element type `T` of a +[`ParticleList`](@ref). + +The mass must be positive: a `mass ≤ 0` throws an `ArgumentError`, and the conversion to `T` +checks again, so a mass that rounds to zero in `T` throws too. The charge may be zero. + +# Example + +```julia +s = Species(:e, -1.0, 0.5) +s.charge # -1.0 +s.mass # 0.5 +``` +""" +struct Species{T} + name::Symbol + charge::T + mass::T + + function Species{T}(name::Symbol, charge, mass) where {T} + m = convert(T, mass) + m ≤ zero(T) && + throw(ArgumentError("the mass of species $name must be positive, got $m")) + new{T}(name, convert(T, charge), m) + end +end + +function Species(name::Symbol, charge, mass) + Species{promote_type(typeof(charge), typeof(mass))}( + name, charge, mass) +end + +Species{T}(s::Species) where {T} = Species{T}(s.name, s.charge, s.mass) + +""" + charge(s::Species) + +Return the charge `q` of the species `s`. +""" +charge(s::Species) = s.charge + +""" + mass(s::Species) + +Return the mass `m` of the species `s`. +""" +mass(s::Species) = s.mass + +""" + species(s::Species) + +Return the species `s` itself, so that `species` reads a species as well as the species of a +particle in a [`ParticleList`](@ref). +""" +species(s::Species) = s diff --git a/test/particle_list.jl b/test/particle_list.jl index 4311cd7..8ade9d8 100644 --- a/test/particle_list.jl +++ b/test/particle_list.jl @@ -101,3 +101,225 @@ plh5 = ParticleList(h5file) @test plh5[1, 1] == pl[1, 1] @test plh5[1, :] == pl[1, :] @test plh5[:, 1] == pl[:, 1] + +# A function barrier with concrete arguments, for the allocation assertions: the first call warms +# the method and the second is measured. +allocations(f::F, a::A) where {F, A} = (f(a); @allocated f(a)) +allocations(f::F, a::A, i::I) where {F, A, I} = (f(a, i); @allocated f(a, i)) + +_allocating(pl, i) = zeros(i) +_chargeweight(pl, i) = chargeweights(pl)[1, i] + +# the species table, the per-particle index, and the accessors +@testset "species table and accessors ($T)" for T in (Float32, Float64) + np = 5 + nd = 3 + vars = (x = 1, v = 2, w = 3) + z = rand(T, nd, np) + specieslist = [Species(:e, -1, T(0.5)), Species(:i, 2, 4)] + speciesindex = [1, 2, 1, 2, 2] + pl = ParticleList(z; variables = vars, species = specieslist, speciesindex = speciesindex) + + # the table is converted to the element type of the list, and the index stays `Int` + @test pl.species isa Vector{Species{T}} + @test pl.speciesindex isa Vector{Int} + @test pl.speciesindex == speciesindex + + # every particle reads its own charge, mass and species + @test [charge(pl, i) for i in 1:np] == T[-1, 2, -1, 2, 2] + @test [mass(pl, i) for i in 1:np] == T[T(0.5), 4, T(0.5), 4, 4] + @test [species(pl, i).name for i in 1:np] == [:e, :i, :e, :i, :i] + @test charge(pl, 1) isa T + @test mass(pl, 1) isa T + @test species(pl, 1) isa Species{T} + + # the accessors infer, and allocate nothing after a warm call + @test @inferred(charge(pl, 1)) == T(-1) + @test @inferred(mass(pl, 1)) == T(0.5) + @test @inferred(species(pl, 1)) == Species(:e, -1, T(0.5)) + @test allocations(charge, pl, 1) == 0 + @test allocations(mass, pl, 1) == 0 + @test allocations(species, pl, 1) == 0 + @test allocations(_allocating, pl, 1) > 0 + + # an `Int32` index works, on the accessors and in the keyword + @test charge(pl, Int32(2)) == charge(pl, 2) + @test mass(pl, Int32(2)) == mass(pl, 2) + @test species(pl, Int32(2)) === species(pl, 2) + @test ParticleList(z; variables = vars, species = specieslist, + speciesindex = Int32[1, 2, 1, 2, 2]) == pl + + # construction throws in three cases + @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist, + speciesindex = [1, 2, 1, 2]) + @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist, + speciesindex = [1, 2, 3, 1, 2]) + @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist) + + # equality compares the species table and the index + @test ParticleList(z; variables = vars, species = specieslist, + speciesindex = speciesindex) == pl + @test ParticleList(z; variables = vars, + species = [Species(:e, -1, T(0.5)), Species(:i, 3, 4)], + speciesindex = speciesindex) != pl + @test ParticleList(z; variables = vars, species = specieslist, + speciesindex = [2, 1, 1, 2, 2]) != pl + + # a list built without a species holds one species, with charge and mass one + pld = ParticleList(z; variables = vars) + @test pld.species == [Species(:default, one(T), one(T))] + @test pld.speciesindex == ones(Int, np) + @test charge(pld, 1) == one(T) + @test mass(pld, 1) == one(T) + @test species(pld, 1).name === :default + @test vec(collect(chargeweights(pld))) == collect(pld.w) + + # an empty species table with np ≥ 1 throws through the out-of-range check + @test_throws ArgumentError ParticleList(z; variables = vars, species = Species[]) + + # an empty list builds, with and without a species + z0 = rand(T, nd, 0) + @test length(ParticleList(z0; variables = vars, species = Species[])) == 0 + @test length(ParticleList(z0; variables = vars)) == 0 + + # the conversion to the list's element type re-checks the mass + if T === Float32 + @test_throws ArgumentError ParticleList(z; variables = vars, + species = [Species(:n, 1.0, 1e-50)], + speciesindex = ones(Int, np)) + end +end + +# the charge-weighted weights, and a deposit that uses them +@testset "chargeweights" begin + np = 5 + x = [0.4, 0.6, 1.4, 1.6, 1.9] + v = zeros(np) + w = [0.25, 0.5, 1.0, 2.0, 4.0] + specieslist = [Species(:e, -1, 0.5), Species(:i, 2, 4)] + speciesindex = [1, 2, 1, 2, 2] + pl = ParticleList(x, v, w; species = specieslist, speciesindex = speciesindex) + + cw = chargeweights(pl) + @test cw isa AbstractMatrix{Float64} + @test size(cw) == (1, np) + # q_i · w_i by hand: (−1)(0.25), (2)(0.5), (−1)(1.0), (2)(2.0), (2)(4.0) + @test collect(cw) == [-0.25 1.0 -1.0 4.0 8.0] + @test @inferred(chargeweights(pl)[1, 1]) == -0.25 + @test cw[1, Int32(2)] == 1.0 + @test allocations(_chargeweight, pl, 1) == 0 + @test_throws BoundsError cw[2, 1] + + # a scalar `w` row and a one-row range agree, on one matrix + z = vcat(x', v', w') + pl_int = ParticleList(z; variables = (x = 1, v = 2, w = 3), species = specieslist, + speciesindex = speciesindex) + pl_rng = ParticleList(z; variables = (x = 1, v = 2, w = 3:3), species = specieslist, + speciesindex = speciesindex) + @test collect(chargeweights(pl_int)) == + collect(chargeweights(pl_rng)) == [-0.25 1.0 -1.0 4.0 8.0] + + # `w` spanning two rows, and a list without `w` + z4 = vcat(x', v', w', w') + @test_throws ArgumentError chargeweights(ParticleList(z4; + variables = (x = 1, v = 2, w = 3:4), species = specieslist, speciesindex = speciesindex)) + @test_throws ArgumentError chargeweights(ParticleList(z; variables = (x = 1, v = 2), + species = specieslist, speciesindex = speciesindex)) + + # the deposit weights each particle by q·w, onto its nearest grid point + grid = [0.5, 1.5] + deposit = zeros(2) + for i in 1:np + j = argmin(abs.(grid .- pl.list[1, i])) + deposit[j] += cw[1, i] + end + # cell 1: (−1)(0.25) + (2)(0.5) = 0.75; cell 2: (−1)(1.0) + (2)(2.0) + (2)(4.0) = 11.0 + @test deposit == [0.75, 11.0] +end + +# the species table and the index survive the HDF5 round trip +@testset "HDF5 round trip of the species table ($T)" for T in (Float32, Float64) + np = 5 + nd = 3 + vars = (x = 1, v = 2, w = 3) + z = rand(T, nd, np) + specieslist = [Species(:e, -1, T(0.5)), Species(:i, 2, 4)] + speciesindex = [1, 2, 1, 2, 2] + pl = ParticleList(z; variables = vars, species = specieslist, speciesindex = speciesindex) + + mktempdir() do dir + f = joinpath(dir, "species.h5") + h5open(f, "w") do file + ParticleMethods.h5save(file, pl) + end + pl2 = ParticleList(f) + @test charge.(Ref(pl2), 1:np) == charge.(Ref(pl), 1:np) + @test mass.(Ref(pl2), 1:np) == mass.(Ref(pl), 1:np) + @test pl2.speciesindex == pl.speciesindex + @test pl2.species == pl.species + end +end + +# a file whose species group is gone loads with the default species +@testset "a file without a species group loads with the default" begin + np = 5 + z = rand(3, np) + vars = (x = 1, v = 2, w = 3) + pl = ParticleList( + z; variables = vars, species = [Species(:e, -1, 0.5), Species(:i, 2, 4)], + speciesindex = [1, 2, 1, 2, 2]) + + mktempdir() do dir + f = joinpath(dir, "stripped.h5") + h5open(f, "w") do file + ParticleMethods.h5save(file, pl) + end + h5open(f, "r+") do file + delete_object(file, "species") + end + pl2 = ParticleList(f) + @test species(pl2, 1).name === :default + @test charge(pl2, 1) == 1.0 + @test mass(pl2, 1) == 1.0 + @test pl2.speciesindex == ones(Int, np) + end +end + +# the species group is written under `path`, beside `particles` +@testset "the species group follows path" begin + np = 5 + z = rand(3, np) + vars = (x = 1, v = 2, w = 3) + plA = ParticleList( + z; variables = vars, species = [Species(:e, -1, 0.5), Species(:i, 2, 4)], + speciesindex = [1, 2, 1, 2, 2]) + plB = ParticleList( + z; variables = vars, species = [Species(:p, 3, 7), Species(:q, -4, 2)], + speciesindex = [2, 1, 2, 1, 2]) + + mktempdir() do dir + f = joinpath(dir, "nested.h5") + h5open(f, "w") do file + ParticleMethods.h5save(file, plA; path = "a") + ParticleMethods.h5save(file, plB; path = "b") + end + a = ParticleList(f, "a") + b = ParticleList(f, "b") + @test charge.(Ref(a), 1:np) == charge.(Ref(plA), 1:np) + @test mass.(Ref(a), 1:np) == mass.(Ref(plA), 1:np) + @test a.speciesindex == plA.speciesindex + @test charge.(Ref(b), 1:np) == charge.(Ref(plB), 1:np) + @test mass.(Ref(b), 1:np) == mass.(Ref(plB), 1:np) + @test b.speciesindex == plB.speciesindex + @test a.species != b.species + + # a nested species group stripped from one path leaves the other alone + h5open(f, "r+") do file + delete_object(file, "a/species") + end + a2 = ParticleList(f, "a") + b2 = ParticleList(f, "b") + @test species(a2, 1).name === :default + @test charge.(Ref(b2), 1:np) == charge.(Ref(plB), 1:np) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index ce6c592..8e894f9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ const GROUPS = isempty(ARGS) ? ["core", "slow"] : ARGS if "core" in GROUPS @safetestset "Aqua" include("quality/aqua.jl") @safetestset "Particles" include("particle.jl") + @safetestset "Species" include("species.jl") @safetestset "Particle Lists" include("particle_list.jl") @safetestset "Defects" include("integration/defects.jl") end diff --git a/test/species.jl b/test/species.jl new file mode 100644 index 0000000..8ed4400 --- /dev/null +++ b/test/species.jl @@ -0,0 +1,39 @@ +using ParticleMethods +using Test + +# a species carries a name, a charge and a mass, and the mass must be positive +@testset "Species rejects a non-positive mass" begin + @test_throws ArgumentError Species(:n, 0.0, 0.0) + @test_throws ArgumentError Species(:n, 1.0, -1.0) + @test_throws ArgumentError Species(:n, 1.0, 0.0) + + s = Species(:n, 0.0, 1.0) + @test s isa Species{Float64} + @test s.name === :n + @test s.charge == 0.0 + @test s.mass == 1.0 +end + +# the constructor promotes the charge and the mass, and the accessors read the fields +@testset "Species promotes its fields" begin + e = Species(:e, -1, 0.5) + @test e isa Species{Float64} + @test charge(e) == -1.0 + @test mass(e) == 0.5 + @test species(e) === e + + i = Species(:i, 2, 4) + @test i isa Species{Int} + @test charge(i) === 2 + @test mass(i) === 4 +end + +# Species{T}(s) converts an existing species, and re-checks the mass in T +@testset "Species converts to another element type" begin + e = Species(:e, -1, 0.5) + @test Species{Float32}(e) isa Species{Float32} + @test charge(Species{Float32}(e)) === -1.0f0 + @test mass(Species{Float32}(e)) === 0.5f0 + + @test_throws ArgumentError Species{Float32}(Species(:n, 1.0, 1e-50)) +end From add4daddc754f4962ceb117cae170ef230eaf900 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sun, 11 Oct 2026 21:08:21 +0200 Subject: [PATCH 2/5] Check the species index below the range, and the dataset types on the round trip --- src/particle_list.jl | 17 ++++++++++------- test/particle_list.jl | 12 +++++++++++- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/particle_list.jl b/src/particle_list.jl index f713226..40a942e 100644 --- a/src/particle_list.jl +++ b/src/particle_list.jl @@ -10,8 +10,10 @@ _sort_ntuple(nt::NamedTuple) = NamedTuple{_sort_names(nt)}(nt) """ ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(), species = nothing, speciesindex = nothing) - ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple()) - ParticleList(x, v, w; parameters = NamedTuple()) + ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple(), + species = nothing, speciesindex = nothing) + ParticleList(x, v, w; parameters = NamedTuple(), species = nothing, + speciesindex = nothing) ParticleList(h5::H5DataStore, path = "/") ParticleList(fpath::AbstractString, path = "/") @@ -26,9 +28,9 @@ field names `list`, `views`, `params`, `particles`, `variables`, `indices`, `spe `species` is the table of particle species, each a [`Species`](@ref) with a charge and a mass, and `speciesindex` gives the species of each particle as an index into that table. A list built -without `species` holds one species with charge and mass one, so `chargeweights(pl)` is `pl.w`. A -list with more than one species needs an explicit `speciesindex` of length `np`, with every index -in `1:length(species)`. +without `species` holds one species with charge and mass one, so `chargeweights(pl)[1, i]` is the +weight `w_i` of particle `i`. A list with more than one species needs an explicit `speciesindex` of +length `np`, with every index in `1:length(species)`. The second form makes a zero list of element type `DT` with `np` particles of `nd` components. The third form stacks positions `x`, velocities `v` and weights `w`, and defines the variables @@ -216,8 +218,9 @@ end chargeweights(pl::ParticleList) Return the charge-weighted weights of the particles in `pl` as a lazy `1 × np` matrix, whose -`[1, i]` element is `charge(pl, i) * pl.w[1, i]`. This is the weight a charge deposit reads, so a -deposit that weights each particle by its charge changes one argument and no signature. +`[1, i]` element is `charge(pl, i) * w_i`, the weight of particle `i`. This is the weight a charge +deposit reads, so a deposit that weights each particle by its charge changes one argument and no +signature. It throws an `ArgumentError` when `pl` has no variable `w`, or when `w` spans more than one row. """ diff --git a/test/particle_list.jl b/test/particle_list.jl index 8ade9d8..7c53518 100644 --- a/test/particle_list.jl +++ b/test/particle_list.jl @@ -149,11 +149,13 @@ _chargeweight(pl, i) = chargeweights(pl)[1, i] @test ParticleList(z; variables = vars, species = specieslist, speciesindex = Int32[1, 2, 1, 2, 2]) == pl - # construction throws in three cases + # construction throws in three cases; the range is checked on both sides @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist, speciesindex = [1, 2, 1, 2]) @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist, speciesindex = [1, 2, 3, 1, 2]) + @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist, + speciesindex = [1, 2, 1, 2, 0]) @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist) # equality compares the species table and the index @@ -257,6 +259,14 @@ end @test mass.(Ref(pl2), 1:np) == mass.(Ref(pl), 1:np) @test pl2.speciesindex == pl.speciesindex @test pl2.species == pl.species + + # the datasets are written in the element type of the list, `index` as `Int` + h5open(f, "r") do file + @test eltype(read(file["species/charge"])) === T + @test eltype(read(file["species/mass"])) === T + @test eltype(read(file["species/index"])) === Int + @test read(file["species/name"]) isa Vector{String} + end end end From 74953f0dd081a34f02eee8d9c371f72914ce34ab Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sun, 11 Oct 2026 22:20:11 +0200 Subject: [PATCH 3/5] Raise the HDF5 floor to 0.16.16 and simplify the species checks h5save writes the species names as a variable-length string dataset, and HDF5 0.16.0 to 0.16.15 cannot read one on Julia 1.10 (Base._memcpy! is undefined there). That broke the Downgrade job, so the floor of HDF5 = "0.16, 0.17" was no longer true. The species index is copied rather than aliased, its length and range are checked in the inner constructor, and a NaN mass is rejected. The default species is a keyword default, the HDF5 reader passes the species keywords only when the group exists, only() replaces _wrow, and the identity method species(::Species) is removed. The HDF5 round-trip tests compare whole lists with ==. Co-Authored-By: Claude Opus 5.5 --- Project.toml | 2 +- src/particle_list.jl | 84 +++++++++++++++++-------------------------- src/species.jl | 15 +++----- test/particle_list.jl | 35 +++++++----------- test/species.jl | 2 +- 5 files changed, 51 insertions(+), 87 deletions(-) diff --git a/Project.toml b/Project.toml index c8eb66d..7413bfd 100644 --- a/Project.toml +++ b/Project.toml @@ -7,5 +7,5 @@ version = "0.1.2" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" [compat] -HDF5 = "0.16, 0.17" +HDF5 = "0.16.16, 0.17" julia = "1.10" diff --git a/src/particle_list.jl b/src/particle_list.jl index 40a942e..303d2bd 100644 --- a/src/particle_list.jl +++ b/src/particle_list.jl @@ -9,10 +9,10 @@ _sort_ntuple(nt::NamedTuple) = NamedTuple{_sort_names(nt)}(nt) """ ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(), - species = nothing, speciesindex = nothing) + species = [Species(:default, 1, 1)], speciesindex = nothing) ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple(), - species = nothing, speciesindex = nothing) - ParticleList(x, v, w; parameters = NamedTuple(), species = nothing, + species = [Species(:default, 1, 1)], speciesindex = nothing) + ParticleList(x, v, w; parameters = NamedTuple(), species = [Species(:default, 1, 1)], speciesindex = nothing) ParticleList(h5::H5DataStore, path = "/") ParticleList(fpath::AbstractString, path = "/") @@ -27,10 +27,11 @@ field names `list`, `views`, `params`, `particles`, `variables`, `indices`, `spe `speciesindex`. `species` is the table of particle species, each a [`Species`](@ref) with a charge and a mass, and -`speciesindex` gives the species of each particle as an index into that table. A list built -without `species` holds one species with charge and mass one, so `chargeweights(pl)[1, i]` is the -weight `w_i` of particle `i`. A list with more than one species needs an explicit `speciesindex` of -length `np`, with every index in `1:length(species)`. +`speciesindex` gives the species of each particle as an index into that table. The table is +converted to the element type of `list`, and the index is copied. A list built without `species` +holds one species with charge and mass one, so `chargeweights(pl)[1, i]` is the weight `w_i` of +particle `i`. A list with more than one species needs an explicit `speciesindex` of length `np`, +with every index in `1:length(species)`. The second form makes a zero list of element type `DT` with `np` particles of `nd` components. The third form stacks positions `x`, velocities `v` and weights `w`, and defines the variables @@ -53,6 +54,10 @@ struct ParticleList{T, ST <: AbstractMatrix{T}, VT <: NamedTuple, PT <: NamedTup @assert isdisjoint(fieldnames(ParticleList), keys(views)) @assert isdisjoint(fieldnames(ParticleList), keys(params)) @assert isdisjoint(keys(views), keys(params)) + length(speciesindex) == size(list, 2) || throw(ArgumentError( + "speciesindex has length $(length(speciesindex)), expected $(size(list, 2))")) + all(in(eachindex(species)), speciesindex) || throw(ArgumentError( + "speciesindex lies outside 1:$(length(species))")) new{eltype(list), typeof(list), typeof(views), typeof(params), typeof(particles), typeof(variables), typeof(indices)}( list, views, params, particles, variables, indices, species, speciesindex) @@ -61,35 +66,24 @@ end function ParticleList( list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(), - species = nothing, speciesindex = nothing) + species = [Species(:default, one(eltype(list)), one(eltype(list)))], + speciesindex = nothing) svariables = _sort_ntuple(variables) views = map(idx_range -> view(list, idx_range, :), svariables) vars = map(v -> eachslice(v; dims = ndims(v)), views) particles = [Particle(p; variables = svariables, parameters = parameters) for p in eachcol(list)] - T = eltype(list) - np = size(list, 2) - stable = species === nothing ? [Species(:default, one(T), one(T))] : - Species{T}.(species) - sind = if speciesindex === nothing - if length(stable) > 1 - throw(ArgumentError( - "a ParticleList with $(length(stable)) species needs an explicit speciesindex")) - end - ones(Int, np) + if speciesindex === nothing + length(species) > 1 && throw(ArgumentError( + "a ParticleList with $(length(species)) species needs an explicit speciesindex")) + index = ones(Int, size(list, 2)) else - convert(Vector{Int}, speciesindex) + index = collect(Int, speciesindex) end - length(sind) == np || throw(ArgumentError( - "speciesindex has length $(length(sind)), expected $np")) - for i in sind - 1 ≤ i ≤ length(stable) || throw(ArgumentError( - "speciesindex contains $i, outside 1:$(length(stable))")) - end - - ParticleList(list, views, parameters, particles, vars, svariables, stable, sind) + ParticleList(list, views, parameters, particles, vars, svariables, + Species{eltype(list)}.(species), index) end function ParticleList(DT::DataType, np::Int, nd::Int; kwargs...) @@ -186,14 +180,14 @@ Return the [`Species`](@ref) of particle `i` in `pl`, that is `pl.species[pl.spe Return the charge of particle `i` in `pl`, the charge of its [`Species`](@ref). """ -@inline charge(pl::ParticleList, i::Integer) = species(pl, i).charge +@inline charge(pl::ParticleList, i::Integer) = charge(species(pl, i)) """ mass(pl::ParticleList, i::Integer) Return the mass of particle `i` in `pl`, the mass of its [`Species`](@ref). """ -@inline mass(pl::ParticleList, i::Integer) = species(pl, i).mass +@inline mass(pl::ParticleList, i::Integer) = mass(species(pl, i)) """ ChargeWeights{T} <: AbstractMatrix{T} @@ -226,15 +220,7 @@ It throws an `ArgumentError` when `pl` has no variable `w`, or when `w` spans mo """ function chargeweights(pl::ParticleList{T}) where {T} haskey(pl.indices, :w) || throw(ArgumentError("the ParticleList has no variable w")) - return ChargeWeights{T, typeof(pl)}(pl, _wrow(pl.indices.w)) -end - -_wrow(w::Integer) = Int(w) - -function _wrow(w) - length(w) == 1 || - throw(ArgumentError("the variable w must span exactly one row, got $(length(w))")) - return Int(first(w)) + return ChargeWeights{T, typeof(pl)}(pl, only(pl.indices.w)) end function ParticleList(h5::H5DataStore, path::AbstractString = "/") @@ -253,19 +239,15 @@ function ParticleList(h5::H5DataStore, path::AbstractString = "/") params = NamedTuple{paramkeys}(paramvals) - species, speciesindex = if haskey(group, "species") - sgroup = group["species"] - names = Symbol.(read(sgroup["name"])) - charges = read(sgroup["charge"]) - masses = read(sgroup["mass"]) - ([Species(n, q, m) for (n, q, m) in zip(names, charges, masses)], - convert(Vector{Int}, read(sgroup["index"]))) + skwargs = if haskey(group, "species") + s = group["species"] + (species = Species.(Symbol.(read(s["name"])), read(s["charge"]), read(s["mass"])), + speciesindex = read(s["index"])) else - (nothing, nothing) + (;) end - ParticleList(particles; variables = variables, parameters = params, - species = species, speciesindex = speciesindex) + ParticleList(particles; variables = variables, parameters = params, skwargs...) end function ParticleList(fpath::AbstractString, path::AbstractString = "/") @@ -290,9 +272,9 @@ function h5save(h5::H5DataStore, p::ParticleList; path::AbstractString = "/") end sgroup = _create_group(group, "species") - sgroup["name"] = String[String(s.name) for s in p.species] - sgroup["charge"] = [s.charge for s in p.species] - sgroup["mass"] = [s.mass for s in p.species] + sgroup["name"] = [String(s.name) for s in p.species] + sgroup["charge"] = charge.(p.species) + sgroup["mass"] = mass.(p.species) sgroup["index"] = p.speciesindex group["particles"] = p.list diff --git a/src/species.jl b/src/species.jl index 0b89399..fc454d2 100644 --- a/src/species.jl +++ b/src/species.jl @@ -7,8 +7,9 @@ A particle species: a name, a charge `q` and a mass `m`, shared by every particl `Species{T}(s::Species)` converts an existing species to the element type `T` of a [`ParticleList`](@ref). -The mass must be positive: a `mass ≤ 0` throws an `ArgumentError`, and the conversion to `T` -checks again, so a mass that rounds to zero in `T` throws too. The charge may be zero. +The mass must be positive: a mass that is not, `NaN` included, throws an `ArgumentError`, and the +conversion to `T` checks again, so a mass that rounds to zero in `T` throws too. The charge may be +zero. # Example @@ -25,7 +26,7 @@ struct Species{T} function Species{T}(name::Symbol, charge, mass) where {T} m = convert(T, mass) - m ≤ zero(T) && + m > zero(T) || throw(ArgumentError("the mass of species $name must be positive, got $m")) new{T}(name, convert(T, charge), m) end @@ -51,11 +52,3 @@ charge(s::Species) = s.charge Return the mass `m` of the species `s`. """ mass(s::Species) = s.mass - -""" - species(s::Species) - -Return the species `s` itself, so that `species` reads a species as well as the species of a -particle in a [`ParticleList`](@ref). -""" -species(s::Species) = s diff --git a/test/particle_list.jl b/test/particle_list.jl index 7c53518..96eed2f 100644 --- a/test/particle_list.jl +++ b/test/particle_list.jl @@ -149,6 +149,12 @@ _chargeweight(pl, i) = chargeweights(pl)[1, i] @test ParticleList(z; variables = vars, species = specieslist, speciesindex = Int32[1, 2, 1, 2, 2]) == pl + # the list copies the index, so a later write to the caller's vector leaves it valid + idx = copy(speciesindex) + pli = ParticleList(z; variables = vars, species = specieslist, speciesindex = idx) + idx[1] = 7 + @test pli.speciesindex == speciesindex + # construction throws in three cases; the range is checked on both sides @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist, speciesindex = [1, 2, 1, 2]) @@ -254,11 +260,7 @@ end h5open(f, "w") do file ParticleMethods.h5save(file, pl) end - pl2 = ParticleList(f) - @test charge.(Ref(pl2), 1:np) == charge.(Ref(pl), 1:np) - @test mass.(Ref(pl2), 1:np) == mass.(Ref(pl), 1:np) - @test pl2.speciesindex == pl.speciesindex - @test pl2.species == pl.species + @test ParticleList(f) == pl # the datasets are written in the element type of the list, `index` as `Int` h5open(f, "r") do file @@ -287,11 +289,7 @@ end h5open(f, "r+") do file delete_object(file, "species") end - pl2 = ParticleList(f) - @test species(pl2, 1).name === :default - @test charge(pl2, 1) == 1.0 - @test mass(pl2, 1) == 1.0 - @test pl2.speciesindex == ones(Int, np) + @test ParticleList(f) == ParticleList(z; variables = vars) end end @@ -313,23 +311,14 @@ end ParticleMethods.h5save(file, plA; path = "a") ParticleMethods.h5save(file, plB; path = "b") end - a = ParticleList(f, "a") - b = ParticleList(f, "b") - @test charge.(Ref(a), 1:np) == charge.(Ref(plA), 1:np) - @test mass.(Ref(a), 1:np) == mass.(Ref(plA), 1:np) - @test a.speciesindex == plA.speciesindex - @test charge.(Ref(b), 1:np) == charge.(Ref(plB), 1:np) - @test mass.(Ref(b), 1:np) == mass.(Ref(plB), 1:np) - @test b.speciesindex == plB.speciesindex - @test a.species != b.species + @test ParticleList(f, "a") == plA + @test ParticleList(f, "b") == plB # a nested species group stripped from one path leaves the other alone h5open(f, "r+") do file delete_object(file, "a/species") end - a2 = ParticleList(f, "a") - b2 = ParticleList(f, "b") - @test species(a2, 1).name === :default - @test charge.(Ref(b2), 1:np) == charge.(Ref(plB), 1:np) + @test ParticleList(f, "a") == ParticleList(z; variables = vars) + @test ParticleList(f, "b") == plB end end diff --git a/test/species.jl b/test/species.jl index 8ed4400..f6651ef 100644 --- a/test/species.jl +++ b/test/species.jl @@ -6,6 +6,7 @@ using Test @test_throws ArgumentError Species(:n, 0.0, 0.0) @test_throws ArgumentError Species(:n, 1.0, -1.0) @test_throws ArgumentError Species(:n, 1.0, 0.0) + @test_throws ArgumentError Species(:n, 1.0, NaN) s = Species(:n, 0.0, 1.0) @test s isa Species{Float64} @@ -20,7 +21,6 @@ end @test e isa Species{Float64} @test charge(e) == -1.0 @test mass(e) == 0.5 - @test species(e) === e i = Species(:i, 2, 4) @test i isa Species{Int} From 2d85f15647170c2efe27b360a26e305dc0eac5da Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sun, 11 Oct 2026 22:21:09 +0200 Subject: [PATCH 4/5] Record the HDF5 floor and the reserved species names in the changelog The eachslice entry returns to its own figure, 184.0 MB, and names the 8 bytes per particle that the species index adds. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e1b99f2..d838cdf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,7 +23,8 @@ first entry is written. asserts, so `test/` was the wrong home: a profiling harness that runs as part of the suite costs a minute and checks nothing. - **`ParticleList` carries a species table and a per-particle species index.** A new `Species` - type holds a name, a charge and a mass. `charge(pl, i)`, `mass(pl, i)` and `species(pl, i)` + type holds a name, a charge and a mass, and rejects a mass that is not positive (`NaN` + included). `charge(pl, i)`, `mass(pl, i)` and `species(pl, i)` read them back for each particle, and `chargeweights(pl)` returns the charge-weighted weights `q_i · w_i` as a lazy `1 × np` matrix, the weight a charge deposit reads. A list built without a species holds one species with charge and mass one, so `q · w` is `w` and every existing @@ -56,8 +57,9 @@ first entry is written. per particle — a redundant representation of the underlying matrix. It now uses `eachslice` over each variable's view, yielding a zero-allocation `AbstractVector` with the same element type and values. Measured at `np = 10⁶` particles, `nd = 6` - components, three variables: allocation reduced from 328.0 MB to 192.0 MB (41.5% - reduction). + components, three variables: allocation reduced from 328.0 MB to 184.0 MB (43.9% + reduction). The species index, which the constructor builds as well, adds 8 bytes per + particle on top of that; see the species entry under New Features. - **`StaticArrays` removed from `[deps]` in Project.toml.** The `using StaticArrays: MVector` import in `src/ParticleMethods.jl` was the only use in `src/` and was removed by the `Particle(DT, len)` change. It remains in `[compat]` and as a test-only dependency in `[extras]` because the test suite @@ -84,6 +86,11 @@ first entry is written. test dependencies are in `test/Project.toml`; `Random` is added, and the tests use `SafeTestsets`. The root `Project.toml` has no `[extras]`, `[targets]` or test-only compat entries. The suite keeps its 81 tests. +- **HDF5 floor raised to 0.16.16** (`[compat]` `"0.16.16, 0.17"`). `h5save` now writes the species + names as a variable-length string dataset (`species/name`). HDF5.jl 0.16.0 to 0.16.15 cannot read + a variable-length string on Julia 1.10, because its loader calls `Base._memcpy!`, which Julia 1.10 + does not define; 0.16.16 guards that call. A Julia 1.10 environment with an older HDF5 resolves to + 0.16.16 or later. ### Bug Fixes @@ -114,3 +121,7 @@ first entry is written. (`AbstractVector`), but this constructor's concrete return type is now correct. This changes equality semantics: `Particle(Float64, 7) == Particle(MVector{7}(zeros(7)))` was `true` before and is now `false`, because `==` requires the state types to match. +- **`species` and `speciesindex` are field names of `ParticleList`, so they are reserved.** The + constructor asserts that variable and parameter names are disjoint from the field names, so a + variable or parameter named `species` or `speciesindex` now fails with an `AssertionError`. Rename + such a variable or parameter. From cce0d1bbf165e52e5271e3e499a5590ce3b1076d Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Sun, 11 Oct 2026 22:30:39 +0200 Subject: [PATCH 5/5] Record the unused chargeweights API as a known issue Co-Authored-By: Claude Opus 5.5 --- KNOWN_ISSUES.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 KNOWN_ISSUES.md diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md new file mode 100644 index 0000000..af238da --- /dev/null +++ b/KNOWN_ISSUES.md @@ -0,0 +1,13 @@ +# Known issues + +### K1 · `ChargeWeights` and `chargeweights` are exported API with no caller. + +- **location:** `src/particle_list.jl:199` +- **evidence:** No code in this package or in its dependents calls `chargeweights`. The type and + the function add about 40 lines of source and about 46 lines of tests. A deposit that weights + each particle by its charge can write `charge(pl, i) * pl.list[wrow, i]` inline: at `np = 10⁶` + the inline form and `chargeweights(pl)[1, i]` measured 692 µs and 695 µs, in one run. Removing + the pair after a release is a breaking change, so the decision is cheapest before the next + release: remove it, or keep it until the first deposit that uses a species charge calls it. +- **kind:** dead code +- **found:** PR #20