Parallelize ClimaCore sparse matrix assembly + polygon intersections#110
Conversation
The standard regridder's parallel COO assembly is now driven by one trait plus two hooks on the intersection operator, each with a default that reproduces today's area-based weight: - `IntersectionReturnStyle(op)`: `OutOfPlaceSingleResult` (default; `op(src, dst) -> area`, the driver stores the triplet) vs `InPlace` (`op(rows, cols, vals, item, src_tree, dst_tree)` pushes its own COO). - `work_items(op, candidate_pairs)`: the units of parallel work (default: one candidate pair each). - `output_matrix_size(op, src_tree, dst_tree)`: the assembled shape (default: dst-cells x src-cells). `intersection_areas` collapses to a single method (`threaded::BoolsAsTypes`) that resolves the return-style trait once and threads it through the shared `_parallel_coo` / `_assemble_chunk` / `_run_and_store!` skeleton; the removed `compute_intersection_areas` becomes the `OutOfPlaceSingleResult` branch. `DefaultIntersectionOperator` uses every default, so the standard path is behaviorally unchanged. The interface is `@public` so downstream packages can define their own parallel-assembling operators. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Both spectral-element assemblers previously built their sparse matrices in a hand-written serial loop. They now plug into the shared parallel `intersection_areas` driver as two `InPlace` operators: - `SEToFVIntersectionOperator`: per (elem, cell) pair, pushes an Nq^2 block of B-weights per intersection polygon (matrix N_fv x N_nodes). - `FVToSEIntersectionOperator`: overrides `work_items` to make the element (not the candidate pair) the unit of work, since the per-element local mass-matrix solve needs all of an element's cells at once (matrix N_nodes x N_fv). Each entrypoint shrinks to constructing an operator and calling `intersection_areas`; the private `_get_candidate_pairs` and the two assembly loops are gone. Assembly at h_elem=32 / 8 threads is ~5.5x faster. Extra kwargs (e.g. `normalize`) are absorbed by the entrypoint signatures and not forwarded, matching the previous behavior. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
On a small cubed sphere, assert the threaded and serial intersection matrices are exactly equal for both SE->FV and FV->SE and both regular and Gilbert (space-filling-curve) element orderings. There are no cross-pair duplicate (row,col) entries in either direction, so chunking cannot change any sum, making this a tight data-race / ordering regression check. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add a "Customizing the weights" section to how_it_works.md (with OutOfPlaceSingleResult and InPlace worked examples), render the seven new public symbols in the API reference, and update the CLAUDE.md architecture notes to describe the trait + two hooks and the ClimaCore operators. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
`GO.area(::Spherical, polygon)` infers `Any` (GeometryOps v0.1.40 boxes a multi-method local closure — JuliaGeo/GeometryOps.jl#407), so spherical_triangle_area returned `Any` and boxed every quadrature weight in accumulate_principled_b (~4000 heap allocs per intersection polygon). Add a `::Float64` barrier on the area call, and build the ring container with a static `SA[...]` to drop its per-triangle allocation. Value-identical (ClimaCore suite green). SE→FV assembly: ~1624 → ~45 allocs/pair (~36x), 380 MiB → 48 MB; FV→SE benefits via the shared accumulator. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Tighten the docstrings introduced in this branch (the intersection-operator interface in intersection_areas.jl and the ClimaCore SE/FV extension) without dropping any content — all PDF equation references and @ref links preserved. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
27f8110 to
9b66683
Compare
|
I've looked over the changes and verified that tests pass locally (including |
|
BTW: I also ran a
|
|
Hmm for that one it's certainly possible that one might try to construct a reusable regridder on two non-overlapping areas (though it would be a bit dumb). The fix should be simple enough that we can just do it there. In general the reason for complexity is twofold:
|
|
Hi @asinghvi17, sorry for the late reply; was at a conference. I can see how very high resolution setups would benefit from this change. If that's a high priority right now, then I'm okay with you merging this in as it doesn't seem to affect our current ClimaCoupler runs. |
|
will fix docs in a follow up PR |
Adds a new interface for the intersection operator, so that ClimaCore only needs to define its intersection operator, not the full sparse matrix assembly.
Also adds some benchmarks for the ClimaCore regridder construction, as well as performance optimisation / improvements there.