Multi-Threading & Garbage Collection in Mera

Threading in Mera: what is automatic, what you control, and how to measure it. Julia 1.10+.

High-performance parallel computing with MERA.jl: leveraging multi-core processors for accelerated astrophysical data analysis

Mera's readers and projection are threaded internally, so most of the benefit needs nothing from you beyond starting Julia with more than one thread. The rest of this page is about the two things that are your decision: how many threads to allow where, and how to thread your own analysis without racing.

1 Turning threading on

Julia starts with one thread. Give it more at startup, then check Mera sees them:

julia -t 8            # explicit count, the safe default
julia -t auto         # every core on the machine

using Mera, Base.Threads
nthreads()            # > 1 means threading is available

That is the whole setup. Mera's readers and projection are threaded internally and need nothing further; the rest of this page is about when that helps and how not to oversubscribe the machine.

#SBATCH --cpus-per-task=16
julia --threads=16,2 --gcthreads=8

echo "Allocated CPUs: $SLURM_CPUS_PER_TASK"
julia -e 'println("Detected CPUs: ", Sys.CPU_THREADS)'

Thread pools and the GC

Julia 1.10+ runs two thread pools and a parallel garbage collector. --threads=8,2 means eight compute threads plus two interactive ones, which keeps the REPL responsive while a long read is running; --gcthreads sets the collector's own threads.

julia --threads=8,2 --gcthreads=4

To see what you ended up with:

using Base.Threads, LinearAlgebra
nthreads(:default)          # compute pool
nthreads(:interactive)      # interactive pool
Threads.ngcthreads()        # garbage collector
BLAS.get_num_threads()      # BLAS has its own pool, see below

BLAS threads multiply with Julia's rather than sharing them, so on a machine where you want at most N busy cores, keep Julia threads × BLAS threads ≤ N:

BLAS.set_num_threads(min(4, nthreads()))

Starting points

These are starting points, not recommendations: measure on your own data, see Measuring.

The three patterns

Everything later on this page is one of these:

2 What Mera threads for you

The threaded functions

Mera Function Threading Architecture:

┌─ MERA FUNCTION CALL ──────────────────────────────────────────┐
│                                                               │
│  gethydro(info; lmax=10, max_threads=4)                      │
│                     ↓                                         │
│  ┌─ PARALLEL FILE LOADING ─┐   ┌─ PARALLEL TABLE CREATION ─┐  │
│  │ Thread 1: amr_001.out01 │   │ Thread 1: :rho column    │  │
│  │ Thread 2: amr_002.out01 │   │ Thread 2: :vx column     │  │
│  │ Thread 3: amr_003.out01 │   │ Thread 3: :vy column     │  │
│  │ Thread 4: amr_004.out01 │   │ Thread 4: :vz column     │  │
│  └─────────────────────────┘   └───────────────────────────┘  │
│                                                               │
│  projection(gas, [:rho, :T, :vx]; max_threads=3)             │
│                     ↓                                         │
│  ┌─ PARALLEL VARIABLE PROCESSING ──────────────────────────┐  │
│  │ Thread 1: Process :rho → density map                   │  │
│  │ Thread 2: Process :T   → temperature map               │  │
│  │ Thread 3: Process :vx  → velocity map                  │  │
│  └─────────────────────────────────────────────────────────┘  │
└───────────────────────────────────────────────────────────────┘

What to expect from more threads

Threads help when there is enough independent work of the right kind. The parallel dimension differs per operation, so the same thread count pays off differently:

The cell-projection case has a practical consequence worth internalising:

# one call, nine variables — the variables run in parallel
projection(gas, [:sd, :T, :vx, :vy, :vz, :σ, :σx, :σy, :σz], :km_s)

# nine calls — each is a single variable, so each runs essentially serially
[projection(gas, v, :km_s) for v in (:sd, :T, :vx, :vy, :vz, :σ, :σx, :σy, :σz)]

Measured on an M2 Pro (12 cores) with mw_L10 output 300, 28.3M cells, the numbers behind this are in Performance:

A single variable stays flat by construction, not through misconfiguration and not because the work is too light. An axis-aligned projection divides its work across the variables you ask for, taking min(max_threads, nthreads(), number_of_variables) threads, so one variable runs on one thread whatever you offer it. The ten-variable case has ten ways to split and gains about 1.6x before saturating.

Off-axis and :exact projections divide by cell instead, so there a single variable does use every thread. See Performance for both measured.

3 How many threads, and max_threads

Threading pays when there is enough independent work of the right kind. It does not pay for a single small calculation, on a memory-starved machine, or when the storage is already the limit.

Do not thread a threaded function without a budget

This is the mistake worth avoiding. Everything in the table above is already threaded, so the moment you put @threads around a Mera call you are nesting two levels of parallelism, and the inner level still defaults to every thread you have.

# 8 outer tasks, each reading with all 8 threads: 64 concurrent readers on one disk
@threads for i in eachindex(snapshots)
    gas = gethydro(info; lmax=10)
    projection(gas, [:rho, :T, :vx, :vy])
end

Julia will not fall over: its threading is composable and it will not spawn 64 OS threads. What it cannot prevent is resource contention, and the resource is almost never the cores:

• memory bandwidth saturates before the cores do on large AMR reads
• network storage is usually faster with fewer concurrent readers, not more
• CPU caches thrash when many memory-heavy tasks interleave

The budget

You have Threads.nthreads() threads. Spend them once:

outer tasks × inner max_threads ≈ nthreads()

@threads runs min(number of items, nthreads()) iterations at a time, so the outer half of that product is something you can work out rather than guess. On eight threads:

# many items: the outer loop owns the threads
@threads for i in eachindex(snapshots)
    gas = gethydro(info; lmax=10, max_threads=1)
    projection(gas, [:rho, :T]; max_threads=1)
end

# one item: no outer loop at all, ask for the variables together
proj = projection(gas, [:sd, :T, :vx, :vy])

For reading many snapshots, use processes rather than threads

The budget above splits threads inside one Julia process. For reading, that is often the wrong axis. Mera's RAMSES reading is allocation bound, and the allocation rate saturates near 1.5 GB/s per process however many threads it is given. Julia's allocator and GC are per-process, so a second process gets its own share while a second thread does not.

Measured on a 32-thread server, reading two snapshots with a 16-thread budget:

1.71x, from halving the threads per read. For comparison, the whole 1 to 16 thread sweep on a single read gained 1.61x. See Performance.

So when you have many snapshots to get through, try one process per snapshot with a modest thread count each, rather than one process working through them with everything.

Two things decide whether that transfers to your setup. Memory: N concurrent reads need N times the peak, so RAM becomes the binding constraint before threads do. And your filesystem: this was measured on local storage, where reading is allocation bound. On a networked filesystem such as Lustre or GPFS, thousands of file opens may make I/O the limit instead, and then several processes contend for one metadata server and can be slower than one. Measure it rather than assume, it takes a few minutes.

Treat it as a starting point, not a law. Reading is I/O bound, so on slow or networked storage fewer concurrent readers often beat the arithmetic: try max_threads=1 with four outer tasks rather than eight. Measure it, see the next section.

4 Measuring

Guessing a thread count is guessing. Two numbers tell you almost everything.

@time gethydro(info; lmax=12)
# 2.345 seconds (1.23 M allocations: 456.7 MiB, 15.2% gc time)

Wall time is what you care about. % gc time is the warning light: above roughly 15% you are allocating too much, and adding threads will make it worse rather than better, because every thread allocates.

To find the right max_threads for one function on your hardware, sweep it:

using BenchmarkTools

for t in (1, 2, 4, 8, Threads.nthreads())
    dt = @belapsed gethydro($info; lmax=12, max_threads=$t)
    println("max_threads=$t → $(round(dt, digits=3)) s")
end

Two things to expect, both normal rather than misconfiguration:

• reading saturates once the storage does, so more threads stop helping and may hurt on network filesystems
• a single axis-aligned projection stays flat, because it parallelises across the variables you ask for, not within one; ask for several in a single call
• off-axis and :exact projections parallelise by cell, so one variable already uses every thread

Published numbers are in Performance.

5 Allocations and the garbage collector

Julia collects garbage while your code runs, and collection pauses every thread. That makes allocation a threading problem: the more you allocate, the less threading buys you. You do not need to know how the collector works, only how to allocate less.

For Mera this is not a general nicety. Reading a snapshot is allocation bound: read time tracks bytes allocated at a steady 1.30 GB/s, measured across a tenfold range of data on a production run, and a full-resolution read allocated 665 GB and spent 62 s collecting. That is why the reading thread sweep flattens where it does, and it is the reason to give the collector its own threads:

julia -t 16                      # 16 compute threads, and 16 GC threads by default
julia -t 16 --gcthreads=8        # 16 compute, 8 for the GC mark phase

GC threads default to the compute thread count, so you already have them. Note that -t N,M does not set them: the second number there is the interactive thread pool. The flag is --gcthreads, or JULIA_NUM_GC_THREADS, and Threads.ngcthreads() reports what you actually got. The reference benchmarks ran 24 compute and 24 GC, the default. It also follows that a recent Julia is worth running: allocator and collector work lands straight on Mera's dominant cost. The package keeps 1.10 as its supported floor so it keeps working, not because it is the version to choose. See Performance for the measurements.

# allocates a temporary array per operation
total = sum(rho .* volume .* factor)

# no temporaries: the arithmetic is fused into the reduction
total = sum(i -> rho[i] * volume[i] * factor, eachindex(rho))

The three habits that matter:

# 1. write into something you allocated once
out = Vector{Float64}(undef, n)
out .= rho .* volume            # the dot fuses and writes in place

# 2. do not grow arrays in a loop
res = Vector{Float64}(undef, n) # not: res = Float64[]; push!(res, x)

# 3. update in place
gas_data .*= 2                  # not: gas_data = gas_data .* 2

If % gc time is still high after that, give the collector its own threads:

julia -t 8 --gcthreads=4

6 Writing your own threaded code

Three patterns cover nearly everything, and they were introduced in section 1.

Outer loop, for many independent items. Thread the loop, and make each call serial so the tasks do not fight over the same disk:

using Base.Threads

outputs = checkoutputs(path, verbose=false).outputs
masses  = Vector{Float64}(undef, length(outputs))

@threads for i in eachindex(outputs)
    info = getinfo(outputs[i], path, verbose=false)
    gas  = gethydro(info; lmax=10, max_threads=1, verbose=false)
    masses[i] = msum(gas, :Msol)
end

Two things to note. @threads for i in eachindex(...), then index inside: @threads needs an indexable range, so for (i, x) in enumerate(xs) throws at run time. And max_threads=1 is there because gethydro is already threaded, see the budget.

Inner kernel, for one dataset and several quantities. Ask for them in one call and Mera threads across them:

proj = projection(gas, [:sd, :T, :vx, :vy])   # not four separate calls

Mixed, when you want to bound both levels:

tasks = [@spawn projection(gas, v; max_threads=2) for v in (:sd, :T)]
maps  = fetch.(tasks)

Keeping it safe

The only rule: never let two threads write the same thing.

# safe: each thread owns one slot
results = Vector{Float64}(undef, n)
@threads for i in 1:n
    results[i] = compute(i)
end

# safe: an atomic counter
total = Atomic{Float64}(0.0)
@threads for i in 1:n
    atomic_add!(total, compute(i))
end

# WRONG: every thread updates the same variable
total = 0.0
@threads for i in 1:n
    total += compute(i)        # races, and the answer changes run to run
end

For anything more structured than a number, use a lock:

using Base.Threads: ReentrantLock, lock
catalogue = Dict{Int,Any}()
lk = ReentrantLock()
@threads for i in eachindex(items)
    r = analyse(items[i])
    lock(lk) do
        catalogue[i] = r
    end
end

7 When threading does not help

Summary

• Start Julia with an explicit thread count. Never -t auto on a shared node.
• Mera's readers and projection are already threaded. You get that for free.
• Thread the outer loop over snapshots, with max_threads=1 inside, or let Mera thread one call over several variables. Not both, uncontrolled.
• Cores are rarely the limit. Storage and memory bandwidth are, which is what max_threads throttles.
• Watch % gc time. Above 15%, allocate less before adding threads.
• Never let two threads write the same thing.