Repository navigation
ParticleList carries a species table, a charge and a mass - #20
Conversation
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.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #20 +/- ##
==========================================
+ Coverage 85.71% 94.69% +8.98%
==========================================
Files 3 4 +1
Lines 98 132 +34
==========================================
+ Hits 84 125 +41
+ Misses 14 7 -7 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
michakraus
left a comment
There was a problem hiding this comment.
Review of add4dad (head of species/part-42)
Verdict: request changes (posted as a comment). One bug: this change makes the HDF5 [compat] floor false, and the Downgrade job shows it. Three correctness risks. And the code can be cut substantially without loss.
CI
- Required checks (classic protection and the ruleset agree: Julia min and Julia 1 on ubuntu, macOS and windows, and Doctests): all 7 green.
- Downgrade - ubuntu-latest (advisory,
continue-on-error): red, and caused by this PR. Julia 1.10.12 with HDF5 v0.16.0 and HDF5_jll v1.10.5+7.Particle Listserrors withUndefVarError: _memcpy! not definedinHDF5._typed_load(::Type{Ptr{UInt8}}, ...)(HDF5 0.16.0src/typeconversions.jl:315). The call comes fromread(sgroup["name"])atsrc/particle_list.jl:258, through the old round trip attest/particle_list.jl:91. The last Downgrade run onmain(run 37980598372, the same Julia and HDF5 v0.16.0) is green, withParticle Listsat 44/44. The newspecies/namedataset is the first variable-length string that the suite reads. HDF5.jl callsBase._memcpy!without a version guard up to v0.16.15, and guards it from v0.16.16 on (@static if VERSION ≤ v"1.10.0-DEV.1390").
Findings
| # | severity | file:line | claim | evidence |
|---|---|---|---|---|
| 1 | bug | Project.toml:10, src/particle_list.jl:258 |
HDF5 = "0.16, 0.17" is now a false floor. On Julia 1.10 with HDF5 0.16.0 to 0.16.15, no file that h5save writes can be loaded, because the read of species/name reaches Base._memcpy!, which Julia 1.10 does not define. Fix: HDF5 = "0.16.16, 0.17", then re-run Downgrade. |
The Downgrade log above. main is green on the same versions. HDF5.jl v0.16.15 src/typeconversions.jl calls Base._memcpy! unguarded; v0.16.16 guards it. |
| 2 | correctness-risk | src/particle_list.jl:82 |
convert(Vector{Int}, speciesindex) returns the vector of the caller itself, so the list aliases it, and any later write to it bypasses the range check. |
With idx = [1, 2, 1, 2, 2] passed as speciesindex: pl.speciesindex === idx is true. After idx[1] = 7, charge(pl, 1) throws BoundsError: attempt to access 2-element Vector{Species{Float64}} at index [7]. Fix: collect(Int, speciesindex). |
| 3 | correctness-risk | src/species.jl:28 |
m ≤ zero(T) passes a NaN mass, though the docstring says that the mass must be positive. |
Species(:n, 1.0, NaN) returns Species{Float64}(:n, 1.0, NaN). Fix: throw unless the mass compares greater than zero, a test that is false for NaN. |
| 4 | correctness-risk | src/particle_list.jl:51-59, 71-90 |
All species validation is in the keyword constructor. The 8-argument inner constructor stores any table and any index without a check. Species keeps its invariant in its inner constructor; ParticleList should do the same, so that no path builds a list whose charge(pl, i) throws later. |
Read: the inner constructor checks only the three name-disjointness assertions. |
| 5 | quality | src/particle_list.jl:198-238, src/ParticleMethods.jl:29, test/particle_list.jl:195-240 |
ChargeWeights, chargeweights and _wrow (41 lines of source, one exported name, 45 lines of tests) are exported API with no caller. The PR body puts the deposit that would use it in a later change. The lazy type is not faster than the product written inline. Remove it, and add it, if it is still wanted, together with the deposit that reads it. The hand-written deposit loop in the test (231-239) tests its own loop, and goes too. |
One process, warm, single thread, Julia 1.13.1, median of 51 calls, machine load about 9: a loop over cw[1, i] against a loop over charge(pl, i) * pl.w[i]: 0.667 against 0.667 µs at np = 10³, 68.96 against 67.67 µs at 10⁵, 695 against 692 µs at 10⁶. |
| 6 | quality | src/particle_list.jl:71-90, 256-268 |
The species block of the keyword constructor (two nothing sentinels, a nested if, a hand loop, the locals stable and sind) and the (nothing, nothing) branch of the reader reduce to keyword defaults and one check in the inner constructor. The sketch below also closes #2 and #4. |
The sketch below, run as a free function, gives the same result as this PR for the default, a missing index with two species, an Int32 index, an index of 0, an empty table with np = 5 and with np = 0, and does not alias the index. |
| 7 | quality | CHANGELOG.md:25-32 |
The entry leaves out the breaking part: species and speciesindex are now reserved, so a list with a variable or a parameter of either name now throws. |
ParticleList(z; variables = (x = 1, v = 2, w = 3), parameters = (species = 1,)) throws AssertionError: isdisjoint(fieldnames(ParticleList), keys(params)) on this head, and builds on main. The same holds for a variable named speciesindex. |
| 8 | quality | test/particle_list.jl:243-335 |
The three HDF5 testsets compare charge, mass, index and table element by element, but this PR extends == to the table and the index. Each round trip is @test ParticleList(f, path) == pl, plus the four element-type checks at 263-268. The charge.(Ref(pl2), 1:np) and mass.(...) lines repeat what pl2.species == pl.species and pl2.speciesindex == pl.speciesindex already imply. |
ParticleList(f) == pl, ParticleList(f, "a") == plA and ParticleList(f, "b") == plB are all true, in Float32 and Float64. |
| 9 | nit | test/particle_list.jl:110, 143 |
The test that allocations(_allocating, pl, 1) is positive tests that zeros(1) allocates. That is Base, not this package. |
Read. |
| 10 | nit | src/species.jl:61, src/particle_list.jl:189, 196 |
species(s::Species) = s is an identity method with no caller. charge(pl, i) and mass(pl, i) read the fields .charge and .mass instead of the accessors charge and mass that this PR defines. |
Read. |
| 11 | nit | src/particle_list.jl:64, 73, 256, 258 |
Local names shadow functions: the keyword species and the reader local species shadow the exported species, and names shadows Base.names. stable reads as an adjective. |
Read. |
| 12 | nit | src/particle_list.jl:261-262, 293-295 |
The reader is Species.(Symbol.(read(s["name"])), read(s["charge"]), read(s["mass"])), and its convert(Vector{Int}, ...) repeats the conversion of the constructor. The writer can use charge.(p.species) and mass.(p.species). |
Read. |
| 13 | nit | CHANGELOG.md:59 |
The eachslice bullet now gives that change 192.0 MB and 41.5 %. That figure contains the 8 MB index of this PR. The eachslice change brought it to 184.0 MB, and the new species bullet already states its own 8 bytes per particle. |
Warm @allocated ParticleList(z; variables) at np = 10⁶, nd = 6, three variables: 184 010 768 B on main, 192 022 688 B on this head. |
| 14 | nit | CHANGELOG.md:52 (adjacent, not from this PR) |
The bullet cites test/particle_list_tests.jl:81, which does not exist; temp.h5 is written at test/particle_list.jl:85. |
Read. |
| 15 | nit | PR body | The body is written for an internal process, not for a reader on GitHub: part numbers (P42, P13, §5V.P42), sections for "the dispatcher" and "advisor decisions", and local script names (mutate.jl, round-1a/probe_blocking.jl). Its "Pre-PR verification" section opens with "Verdict: stop and ask — 2 unresolved verification gates". It names HDF5 0.16.0 as the floor that was checked for delete_object, but no run at that floor read the new string dataset (#1). Its "cold" allocation figures are warm figures: a cold first call measures 279.8 MB on main and 290.5 MB here, compilation included. The delta it gives is correct. |
The PR body; the measurements under "Verified". |
The simplification for #2, #4 and #6
# inner constructor: the one place that every path goes through
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))
sind = collect(Int, speciesindex)
length(sind) == size(list, 2) && all(in(eachindex(species)), sind) ||
throw(ArgumentError("speciesindex needs one entry per particle, each an index into species"))
new{...}(list, views, params, particles, variables, indices,
Species{eltype(list)}.(species), sind)
endfunction ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(),
species = [Species(:default, one(eltype(list)), one(eltype(list)))],
speciesindex = length(species) ≥ 2 ?
throw(ArgumentError("more than one species needs a speciesindex")) :
ones(Int, size(list, 2)))
# ... views, vars, particles as before ...
ParticleList(list, views, parameters, particles, vars, svariables, species, speciesindex)
end# reader
kw = 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, parameters = params, kw...)This removes the nothing sentinels, the nested if, the hand loop and the (nothing, nothing) branch. It copies the index and checks it in one place.
Verified
- Downgrade attribution: the job log on this head, against the Downgrade job of run 37980598372 on
main(Julia 1.10.12, HDF5 v0.16.0, green). HDF5.jl v0.16.0 and v0.16.15 callBase._memcpy!unguarded; v0.16.16 and v0.17.0 guard it. - Required checks:
branches/main/protectionandrules/branches/mainlist the same 7 contexts, and all 7 pass on this head. Downgrade is not required. - Construction allocation: the rise on this head is exactly the index and the table. Warm, at np = 10, 10³, 10⁵ and 10⁶, head minus
mainequals@allocated ones(Int, np)plus 80 B (224 B = 144 + 80 at np = 10; 8 011 920 B = 8 011 840 + 80 at 10⁶). ParticleList(f, path) == plafterh5save, at the root and under two nested paths, in Float32 and Float64.- JuliaFormatter (
sciml) reports all six changed.jlfiles as formatted. - No closing issue references.
- The local checkouts of VlasovMethods and ReducedBasisMethods use no variable or parameter named
speciesorspeciesindex.
Not checked
- A local run of HDF5 0.16.16 on Julia 1.10: the probe environment could not download artifacts. The bump in #1 needs the Downgrade job to confirm it.
- The full test suite was not run locally; CI runs it.
- A name clash of the new exports
charge,massandspecieswith other packages that a user loads together with this one.
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 <[email protected]>
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 <[email protected]>
Review findings: resolutionThis resolves every finding of the review on
|
Co-Authored-By: Claude Opus 5.5 <[email protected]>
ParticleList carries a species table, a charge and a mass
First-class charge, mass and species for the particle lists: a
Speciestype, a species table and aper-particle index on
ParticleList, accessors that infer and allocate nothing, a lazycharge-weighted weight matrix, and the HDF5 round trip of the new data.
What this does
Species{T}(newsrc/species.jl)name::Symbol,charge::T,mass::T;Species(name, charge, mass)promotes the twonumbers to their common type.
NaNincluded, throwsArgumentError; a zero charge is allowed.The check lives in the inner constructor, so
Species{T}(s::Species)and the conversion of aFloat64table into aFloat32list check again, and a mass that rounds to zero inTthrows.charge(s::Species)andmass(s::Species), both exported.ParticleList(src/particle_list.jl)species::Vector{Species{T}}withT = eltype(list), andspeciesindex::Vector{Int}of length
np. One matrix, one list: the splitting method copiesparticles.zas one state andParticleDistributionwraps oneParticleList, so a tuple of lists, one per species, would breakboth.
speciesandspeciesindexkeywords reach the(DT, np, nd)and(x, v, w)forms throughkwargs, so no caller's shape changes and VlasovMethods needs no edit. The keyword constructorcopies
speciesindex, so a later change to the caller's vector does not reach the list.species(pl, i),charge(pl, i)andmass(pl, i):@inline, inferred, and no allocationafter a warm call. They take any
Integerindex.Base.:(==)compares the two new fields aswell.
ArgumentErrorwhenlength(speciesindex) ≠ npor when an indexlies outside
1:length(species), so the positional form checks too. The keyword form also throwsfor more than one species without a
speciesindex.Species(:default, one(T), one(T))withspeciesindex = ones(Int, np), soq · w == wbit for bit and every existing result isunchanged. The index costs 8 bytes per particle at construction.
speciesandspeciesindexare now field names, so a variable or parameter with either namefails the constructor's
AssertionError. The CHANGELOG lists this under Breaking Changes.chargeweights(pl)(exported;ChargeWeightsis not)ChargeWeights{T} <: AbstractMatrix{T}of size1 × npwhose[1, i]element ischarge(pl, i) * w_i. That is the shape and the value a merged_deposit!(potential, vals, x, w)reads as
w[1, i], so a later switch of a deposit toq · wchanges one argument at each callsite, needs no buffer and leaves the signature alone. It reads
pl.list[r, i]at the one row ofw, so a scalarwrow and a one-row range agree.ArgumentErrorwhen the list has no variablew, or whenwspans more than one row.HDF5
h5savewrites a groupspeciesunderpath, besideparticles, withname(
Vector{String}),chargeandmass(inT) andindex(Vector{Int}); the reader reads itfrom the same group.
HDF5 = "0.16.16, 0.17"). The species names are the firstvariable-length strings the package reads, and HDF5.jl 0.16.0 to 0.16.15 cannot read one on
Julia 1.10: their loader calls
Base._memcpy!, which Julia 1.10 does not define. The Downgradejob fails on 0.16.0 and passes on 0.16.16.
Scope
This PR does not edit VlasovMethods. The switch of a deposit to a species weight belongs to the
later change that first deposits with a species charge.
Tests
test/species.jl(new): the non-positive andNaNmass, the promotion ofchargeandmass,the conversion to another element type and its check.
test/particle_list.jl:species table and accessors (Float32)and(Float64),chargeweights,HDF5 round trip of the species table (Float32)and(Float64),a file without a species group loads with the default,the species group follows path. They cover anInt32index, an indexbelow and above the range, a copied index, a mass that rounds to zero in
Float32, the scalar andthe range form of
w, a two-roww, and the element type of the datasets on disk.Known issues added
ChargeWeightsandchargeweightsare exported API with no caller. Remove the pairbefore the next release, or keep it until the first deposit that uses a species charge calls it.
🤖 Generated with Claude Code