Materials: a registry and two distributions, and a JIT constant-collision fix

[lmoresi]

Started as two particle demonstrations; ended up as a material system and a silent correctness bug in the JIT.

The bug first

Now split out as #717, a bugfix/ branch off development, since it is a correctness fix independent of everything else here. The commit below (d53d0038) is the same content; the branches merge cleanly except for tests/test_0103_jit_rampable_constants.py, where this branch's end-to-end test uses the MaterialSwarm API introduced here — take this branch's version of that file.

A layered viscosity built from two ViscousFlowModels returned the uniform-viscosity answer: L2 2.819e-1 against the exact layered Couette profile, and 1.502e-10 against plain linear shear. The composed flux was correct and so was the constants manifest — (0, \eta = 1.0), (1, \eta = 1000.0) — and only the emitted C was wrong.

_JITConstant is a plain sympy.Symbol subclass that never adopted the disambiguation UWexpression uses (docs/developer/design/SYMBOL_DISAMBIGUATION_2025-12.md), and it had neither half of it. It constructed through Symbol.__new__, which is cached by name:

a, b = _JITConstant(0, name='same'), _JITConstant(1, name='same')
a is b                    # True
a._ccodestr, b._ccodestr  # constants[1], constants[1]

The second placeholder was the first object, and assigning its _ccodestr overwrote the first one's. Every occurrence then rendered as one constants[] slot, so the blend collapsed to (phi_0 + phi_1) * constants[k]. Every ViscousFlowModel calls its viscosity \eta, so any model with two of them was affected, not just multi-material.

Fixed structurally: construct via Symbol.__xnew__ to bypass the cache, and put the slot index in _hashable_content, so identity is the slot rather than the name. The same solve now gives 1.807e-7. Separately, the manifest sort tie-broke on str(expr) — a UWexpression's current value — so two same-named constants could swap slots when a parameter changed; it now tie-breaks on instance_number.

Regression tests at both levels in test_0103_jit_rampable_constants.py.

Materials

A model script should not write a level set, a mask, or a blend.

materials = uw.swarm.MaterialSwarm(mesh, fill_param=3)
materials.add("mantle", shear_viscosity_0=1.0,   density=3300)
materials.add("slab",   shear_viscosity_0=1.0e3, density=3400)
materials["slab"] = mesh.X[1] > 0.53

stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.materials = materials
stokes.bodyforce = -materials.density * mesh.CoordinateSystem.unit_e_1

solver.materials sets every constitutive-model parameter the materials declare and the model recognises, by name. One the model does not own (density) is not pushed anywhere; it is a blended symbol for the script to use where it belongs.

What a material is and where it is are kept apart — which is what the existing materials register was reaching for. MaterialRegistry holds the definitions: mesh-free, shareable, exportable, and its values may be numbers, quantities, or laws (eta_0 * sympy.exp(-T.sym[0])). Two distributions consume one registry:

	MaterialSwarm	MaterialRegions
material moves	yes	no
level sets	sampled from particles	exact, 0 or 1
population control	needed	none
cost per step	a proxy fill	nothing, built once

MaterialRegions takes a gmsh physical group, a mesh label, a (label, value) pair, a geometric condition resolved at the integration points, a boolean array, or a callable.

The partition of unity earns its place precisely when a property is a law: there is then no number to store at an integration point, and the only way to combine N expressions into one symbol the assembler can compile is the weighted sum.

MaterialRegistry keeps create_material / set_property / get_property / export_config / import_config / add_callback and the create_standard_* helpers. Removed: assign_to_region and evaluate_property_field, which returned numpy from region IDs and could never reach a weak form — MaterialRegions is that idea, finished. IndexSwarmVariable stays as the machinery and stays importable, but is no longer the documented route.

A declared property that no model recognises and nothing reads is reported at solve time only when it is a close match to a parameter the model does have: viscocity is flagged against viscosity, density is silent.

proxy_sampling

proxy_location says where a swarm variable's proxy lives; the new proxy_sampling says what each point reads.

• "reconstruct" — a weighted fit over the nnn nearest particles, linear-exact. What a smooth field wants (the default).
• "share" — the mean over the particles this point speaks for: those whose nearest integration point within their own cell is this one. The cell-restricted Voronoi share.

The share is what a history wants — stress, damage, accumulated strain — because it is bounded by the particle values, so it cannot invent a stress the swarm never held, while the reconstruction overshoots a discontinuity and gets worse as particles are added:

particles/cell	reconstruct range	share range
3	−0.051 … +1.137	0.000 … 1.000
8	−0.068 … +1.108	0.000 … 1.000
15	−0.109 … +1.097	0.000 … 1.000

There is deliberately no "nearest" for a plain SwarmVariable: sampling one particle's value whole is the material mapping, and putting the property field on the particles hands the solver an answer where it needs a constitutive law. Asking for it raises, and the message says so.

IndexSwarmVariable now defaults to proxy_location="integration_points". Nodes measured worst everywhere (layered Couette 8.0e-2 against 1.8e-7) and its smear is about one cell wide however many particles are added — a property of the basis, not of the swarm.

At a cut-cell interface, nearest sampling converges to a floor set by the rule at ~8 particles per cell (3.48e-2 → 1.85e-2 → 1.85e-2 at fill 3/8/15); past that, refine the mesh. The share gives fractional masks there, and then the mixing rule matters: createMask is arithmetic (Voigt) and does not converge (3.53e-2 flat), while the same masks blended harmonically (Reuss) do (1.83e-2 → 1.24e-2). materials.mixing(shear_viscosity_0="harmonic") selects it.

Locating particles is the expensive half of any cell-local operation (18.8 ms for 32,912 particles against 3.2 ms for the share itself), so Swarm._owning_cells() caches it and repopulate's census reads the same cache. Lagrangian_Swarm forwards proxy_sampling, so a viscoelastic stress history can ride on the share.

Repopulation keeps a label a label

Swarm.repopulate gave new particles a reconstruction of their neighbours' values, which is wrong for an integer material index: the average of two labels is not a label. Integer variables now take their nearest neighbour's value whole; nearest= names any other variable to treat the same way.

Demonstrations and docs

docs/advanced/particle-population-and-materials.md with figures, and runnable scripts in docs/examples/utilities/intermediate/: population control in an extending box, and the layered material problem. docs/api/materials.md and docs/developer/subsystems/integration-point-variables.md updated.

Tests

test_0071 (16), test_0072_material_swarm.py (13), test_0073_material_regions.py (7), test_0103 (+2). Full level_1 and tier_a suite: 1223 passed, 3 skipped, 1 xfailed.

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G

[Copilot]

🟡 Changes recommended

There are a few correctness/documentation issues in the new code paths (starved-rank guard for integration-point mapping, nearest= handling for string inputs, and broken/incorrect doc references) that should be fixed before approval.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR extends UW3’s particle material workflow by (1) letting IndexSwarmVariable project material masks at integration points or as per-cell fractions, and (2) ensuring Swarm.repopulate() preserves integer “labels” via nearest-neighbour copying. It also adds two runnable demonstrations and a focused test suite to validate boundedness, partition-of-unity, and accuracy for a layered-viscosity benchmark.

Changes:

• Add proxy_location={"nodes","integration_points","cells"} to IndexSwarmVariable, including direct particle→integration-point mapping and a per-cell polynomial fraction path.
• Add nearest= support to Swarm.repopulate() and automatically treat integer-valued swarm variables as “nearest-copy” during repopulation.
• Add advanced docs, two intermediate examples, and a new level-1/tier-a test module covering the new behaviors.

File summaries

File	Description
src/underworld3/swarm.py	Implements new IndexSwarmVariable.proxy_location behaviors and label-preserving repopulation logic.
tests/test_0071_index_swarm_proxy_location.py	Adds tests for boundedness/partition-of-unity, integration-point exactness, layered Couette accuracy, repopulation label preservation, and population control behavior.
docs/examples/utilities/intermediate/Ex_Swarm_Population_Control.py	Demonstrates population control in a pure-shear extending box.
docs/examples/utilities/intermediate/Ex_Swarm_Material_Index.py	Demonstrates layered-viscosity accuracy improvements from integration-point material mapping.
docs/advanced/particle-population-and-materials.md	New advanced documentation page explaining population control and material mapping, with figures and links to examples.
docs/advanced/index.md	Links the new advanced documentation page from the advanced docs index.

Review details

Suppressed comments (1)

docs/advanced/particle-population-and-materials.md:116

• The documentation points to docs/examples/utilities/Ex_Swarm_Material_Index.py, but the example added in this PR lives under docs/examples/utilities/intermediate/. This broken path will confuse readers and makes the doc page harder to follow.

`docs/examples/utilities/Ex_Swarm_Material_Index.py`. Two viscosity layers,
1 and 1000, carried as a material index and driven from the top. With the
interface on mesh edges the exact velocity is piecewise linear and lies in the
P2 velocity space, so the only error in the solve is how the material is

• Files reviewed: 6/9 changed files
• Comments generated: 4
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[Copilot]

🟡 Changes recommended

proxy_sampling="share" is currently blocked by the local_size <= 1 guard in _rbf_to_meshVar, which can leave stale integration-point proxies on ranks holding exactly one particle.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Review details

• Files reviewed: 22/26 changed files
• Comments generated: 1
• Review effort level: Lite

[lmoresi]

Review from the 2026-09 style/clutter audit. Findings against the Style Charter, scoped to lines this PR adds.

1. swarm.py:2953 — a self-assignment whose comment claims something it does not do.

for var in self._meshLevelSetVars:
    var.data[:, 0] = var.data[:, 0]          # keep, but write collectively

Reading a local array and writing the same values back is not a collective. If
the starved-rank path genuinely needs one — and this is the empty-rank class we
have been chasing for months — then this line is silently skipping it, and the
level-set variables on a starved rank are left in whatever state the previous
step put them in while the other ranks move on. If no collective is needed, the
statement is dead and the comment is wrong. Either way it cannot stay as
written (§4: a comment states the constraint the code cannot show; §5: no dead
code).

Worth saying which it is in the PR: the answer decides whether this is a tidy-up
or a correctness fix.

2. swarm.py:2947 — uw.mpi.rank in a user-visible warning.

f"IndexSwarmVariable proxy update: rank {uw.mpi.rank} holds "

§11: the user should never see parallel/MPI calls; there is a wrapper and it
works better. uw.pprint / uw.selective_ranks exist for exactly this.

3. Flat .data where §7 mandates .array — materials.py:609, 611, 644;
swarm.py:2978, 2987, 2996; swarm_materials.py:214;
_function.pyx:851. These compute or write external values, so they are not
the one sanctioned flat-.data use (a raw variable-to-variable copy inside the
non-dimensionalisation boundary). Filed as #722 with the merged instances,
since the pattern is wider than this PR.

4. The tests this PR adds sit in the test_006* band, which no glob in
scripts/test.sh reaches — see #721. They have not run in CI.

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review round. Nine confirmed defects plus two ergonomic traps, all fixed with regression tests that reproduce them against the previous revision.

The worst one killed the job. getStratumIS() for a value outside a label's live set hard-aborts PETSc — no exception, no traceback, every rank gone — and the guard if stratum is not None never fired, because petsc4py returns a live IS wrapping a NULL handle. The API's own error message steered users straight in ("say which with materials[name] = (label, value)"), and on a UW3 box the only valid value is 99999, so every guess crashed. The repo documents this hazard twice, in its own words; I had not looked.

Label values are rank-local, which made a two-valued label fail loudly in serial and silently paint the union of both zones at np=2. The value set is now reduced across ranks before anything branches on it — identical behaviour at np=1/2/3.

The rest: two distributions sharing a name silently shared level sets (detected by counting what the mesh gained, not by predicting sanitised names); delete_material/registry.add after a build re-pointed the blend; hasattr raised KeyError and triggered the collective build; a property sharing a name with a real attribute was silently unreachable; m["x"] = A; m["x"] = B gave the union; harmonic mixing divided by a zero value; mixing() accepted undeclared names; units were resolved at add() so the same declaration meant numbers 28 orders apart; solver.materials = None kept pushing.

Style charter: model renamed (§3), the library source moved to .array (§7), and both examples rebuilt on uw.Params/uw.pprint with public API and global diagnostics (§11/§1) — they now run correctly on more than one rank.

The population-control example changed substance. Its claim rested on a scatter plot; three field-level diagnostics I tried first did not discriminate at all, because the nearest-particle mapping degrades gracefully. Starvation wrecks the per-cell fit:

proxy_location	population control	layer area vs exact
cells	on	1.03x
cells	off	5.03x
integration_points	on	1.06x
integration_points	off	1.07x

So the honest claim is narrower: population control is cheap and always safe, but "cells" is where it is necessary.

Full level_1 and tier_a: 1230 passed. Material tests 43/43 at np=1 and np=2.

Two things left for a maintainer decision, both noted in the thread: the .data/.array conflict between the charter, CLAUDE.md and the CI gate's own remediation string; and whether the new tests should enter as tier_b rather than tier_a.

Underworld development team with AI support from Claude Code