ParticleList carries a species table, a charge and a mass

CHANGELOG.md

@@ -22,6 +22,15 @@ 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, 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
+  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
 
@@ -49,7 +58,8 @@ first entry is written.
   `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%
-  reduction).
+  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
@@ -76,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
 
@@ -106,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.

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

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"

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

src/particle_list.jl

@@ -8,9 +8,12 @@ _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(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple())
-    ParticleList(x, v, w; parameters = NamedTuple())
+    ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(),
+                 species = [Species(:default, 1, 1)], speciesindex = nothing)
+    ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple(),
+                 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 = "/")
 
@@ -20,7 +23,15 @@ 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. 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
@@ -35,24 +46,44 @@ 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))
+        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)
+            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 = [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)]
-    ParticleList(list, views, parameters, particles, vars, svariables)
+
+    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
+        index = collect(Int, speciesindex)
+    end
+
+    ParticleList(list, views, parameters, particles, vars, svariables,
+        Species{eltype(list)}.(species), index)
 end
 
 function ParticleList(DT::DataType, np::Int, nd::Int; kwargs...)
@@ -84,7 +115,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 +168,61 @@ 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) = 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) = mass(species(pl, i))
+
+"""
+    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) * 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.
+"""
+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, only(pl.indices.w))
+end
+
 function ParticleList(h5::H5DataStore, path::AbstractString = "/")
     group = h5[path]
 
@@ -151,7 +239,15 @@ function ParticleList(h5::H5DataStore, path::AbstractString = "/")
 
     params = NamedTuple{paramkeys}(paramvals)
 
-    ParticleList(particles; variables = variables, parameters = params)
+    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
+        (;)
+    end
+
+    ParticleList(particles; variables = variables, parameters = params, skwargs...)
 end
 
 function ParticleList(fpath::AbstractString, path::AbstractString = "/")
@@ -175,6 +271,12 @@ function h5save(h5::H5DataStore, p::ParticleList; path::AbstractString = "/")
         params[string(key)] = p.params[key]
     end
 
+    sgroup = _create_group(group, "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
 end
 

src/species.jl

@@ -0,0 +1,54 @@
+"""
+    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 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
+
+```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