Particle transport: repopulation, mid-time swarm step, viscoelastic stress history on particles; conditioning guard in the cells fit - #713
Conversation
… ill-conditioned cells Swarm.repopulate(min_per_cell, max_per_cell, values, order) takes the per-cell census and refills a starved cell from its own lattice at the points farthest from the particles present. A new particle takes the bounded Shepard reconstruction of every swarm variable from its nearest neighbours (order=1 for linear-exact; measured extrapolating to 100 on a field bounded by 1 in a starved corner), or a supplied value. swarm.population_control = dict(...) runs it at the end of every advection(), before the next fit. Collective (the domain test reduces). The cells fit routes a cell whose Gram matrix condition number exceeds cond_max (1e6) to the patch fit whatever its count, and a flat patch keeps its mean: particles the advection clamps onto a wall slide along it as a line, the P2 fit of a line is singular (condition 1e300 at 92 particles), and that garbage grew by 1e12 in ten steps through the PIC read-back. Rotating Gaussian, PIC, ten particles per cell, C = 0.25: population control takes the L2 error from 1.6e-2 to 8.8e-3, level with the integration-point history; the untapered box with wall in- and outflow gives the same answer clamped or with exiting particles deleted. The cap (max_per_cell, closest-pair removal) costs accuracy and is off by default. Tests: test_0068 (refill exact for linear fields with order=1, value override and cap, population control under rotation), a collinear-cell test in test_0067; the history-before-advection test made rank-safe. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
There was a problem hiding this comment.
🟡 Changes recommended
Swarm.repopulate() adds particles via DMSwarm.addNPoints() but misses the PETSc<3.24 empty-local-swarm off-by-one workaround already used in populate(), which can break repopulation on ranks with zero local particles.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Pull request overview
This PR hardens particle-in-cell workflows that use the cells-based swarm proxy by (1) repopulating starved cells after advection and (2) guarding the per-cell polynomial fit against ill-conditioned (e.g., collinear) particle configurations, preventing singular fits and blow-ups.
Changes:
- Added
Swarm.repopulate()plusswarm.population_controlto automatically enforce per-cell particle counts afteradvection(). - Added a condition-number guard to
CellPolynomialProjector.fit()to route ill-conditioned cells to the patch fit (and flatten the patch fit to a mean when needed). - Added/updated tests and extended developer documentation for repopulation and conditioning behavior.
File summaries
| File | Description |
|---|---|
| tests/test_0068_swarm_repopulation.py | New level_1/tier_a tests covering repopulation behavior, value reconstruction/overrides, and population control under rotation. |
| tests/test_0067_integration_point_proxy.py | Makes history-before-advection test rank-safe and adds a regression test for collinear-particle conditioning in the cells proxy. |
| src/underworld3/utilities/cell_polynomial_projection.py | Adds cond_max conditioning guard for dense-cell fits and a “flat patch → mean only” fallback. |
| src/underworld3/swarm.py | Introduces Swarm.repopulate() and population_control hook executed at the end of advection(). |
| docs/developer/subsystems/integration-point-variables.md | Documents repopulation/population_control usage and the conditioning guard rationale/behavior. |
Review details
- Files reviewed: 5/5 changed files
- Comments generated: 2
- Review effort level: Lite
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| d = np.linalg.norm(P[:, None, :] - P[None, :, :], axis=2) | ||
| np.fill_diagonal(d, np.inf) | ||
| nearest = d.min(axis=1) | ||
| surplus = int(npc[c] - max_per_cell) |
There was a problem hiding this comment.
Fixed in the follow-up commit: the cap path takes each particle's nearest-neighbour distance from one kd-tree query (k=2) instead of a pairwise matrix per cell, and repopulate carries the same PETSc<3.24 first-add allocation workaround as populate.
| self.dm.finalizeFieldRegister() | ||
| self.dm.addNPoints(n_new) | ||
| coords = self.dm.getField("DMSwarmPIC_coor").reshape((-1, dim)) |
There was a problem hiding this comment.
Fixed in the follow-up commit: the cap path takes each particle's nearest-neighbour distance from one kd-tree query (k=2) instead of a pairwise matrix per cell, and repopulate carries the same PETSc<3.24 first-add allocation workaround as populate.
swarm.advection(..., midtime_velocity=True) takes the mid-point velocity
at the mid time, 1.5 v^n - 0.5 v^{n-1}, from a CharacteristicTrace the
swarm owns (previous velocity cached by evaluation at the nodes at the end
of each call); characteristics=<trace> uses a solver's shared trace. On a
rotation whose rate ramps linearly, ten frozen-velocity steps miss 0.05
rad, the mid-time steps under 0.008. Off by default (a steady flow is
unchanged), ignored under substepping. The per-step advection prints are
now gated on swarm.verbose.
Underworld development team with AI support from Claude Code
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
…kes solver's DFDt Stokes(..., DFDt=Lagrangian_Swarm(swarm, zeros, SYM_TENSOR, degree=1, continuous=False, proxy_location="cells", step_averaging=1)) carries the stress history on the particles and reads it at the integration points through the cells proxy; after each solve the new stress is evaluated at the particles (the Maxwell update is a local ODE: no projection back, no null space). Plumbing: the constitutive-model setter takes the viscoelastic order from a supplied history (the solver's own _order is set only when it builds one); the solve leaves the nodal projection and shift to the nodal history (isinstance guard) and lets the swarm manager update itself; Lagrangian_Swarm.update_pre/post_solve accept the nodal manager's extra keywords, and update_post_solve evaluates every component of the new field at the particles BEFORE shifting its chain, since a stress expression reads the history it is about to overwrite (audit SWARM-06). Maxwell shear box (test_0070): order 1 within 5% after 20 steps at dt = 0.1 t_r, order 2 within 1%, particle and nodal histories agree to 0.2%; two ranks green; the nodal VEP regressions (test_1052, test_1051) unchanged. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
…es, recorded as open Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
…; PETSc<3.24 allocation workaround in repopulate Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
Follows #707. The particle-transport line on the cells proxy: repopulation (and the failure it uncovered), the swarm step at the mid time, and the viscoelastic stress history carried on particles.
Swarm.repopulate and population control
Swarm.repopulate(min_per_cell=None, max_per_cell=None, values=None, nnn=None, order=0): per-cell census from the strict locator; a starved cell is refilled from its own lattice (the pointspopulateuses) at the points farthest from the particles present. New particles take the bounded Shepard reconstruction of every swarm variable from their nearest neighbours (order=1for linear-exact, which extrapolates in a starved corner), or a supplied value (constant or callable, an inflow datum). Collective.swarm.population_control = dict(...)runs it at the end of everyadvection(), before the next fit.max_per_cellthins over-full cells (closest-pair removal); measured it costs accuracy and is off by default.Conditioning guard in the cells fit
Particles the advection clamps back onto a wall slide along it as a line, and the P2 fit of a collinear set is singular whatever its count (condition number 1e300 at 92 particles in one wall cell); the garbage grew by 1e12 in ten steps through the PIC read-back.
CellPolynomialProjector.fitnow routes a cell whose Gram matrix condition number exceedscond_max(1e6) to the patch fit, and a flat patch keeps only its mean.Swarm step at the mid time
swarm.advection(..., midtime_velocity=True)evaluates the RK2 mid-point velocity at the mid time, 1.5 v^n - 0.5 v^{n-1}, from aCharacteristicTracethe swarm owns, or from a solver's shared trace (characteristics=). Ten steps on a linearly ramping rotation: frozen velocity misses 0.05 rad, mid-time under 0.008. Off by default; a steady flow is unchanged. The per-step advection prints are gated onswarm.verbose.Viscoelastic stress history on particles
Stokes(..., DFDt=Lagrangian_Swarm(swarm, zeros, SYM_TENSOR, degree=1, continuous=False, proxy_location="cells", step_averaging=1))carries the stress history on the particles and reads it at the integration points through the cells proxy; after each solve the new stress is evaluated at the particles (the Maxwell update is a local ODE, so there is no projection back and no null space). Plumbing: the constitutive-model setter takes the viscoelastic order from a supplied history; the solve leaves the nodal projection and shift to the nodal history;Lagrangian_Swarmaccepts the nodal manager's keywords and evaluates every component of the new field before shifting its chain (the stress expression reads the history it overwrites). Maxwell shear box: order 1 within 5%, order 2 within 1%, particle and nodal histories agree to 0.2%.Measured (rotating Gaussian, PIC, ten particles per cell, C = 0.25, one revolution): population control takes the L2 error from 1.6e-2 to 8.8e-3, level with the integration-point history at 9.1e-3. The untapered box, where every wall has an inflow and an outflow segment, gives the same answer whether exiting particles are clamped or deleted (
mesh.return_coords_to_bounds = None).Tests:
test_0068_swarm_repopulation.py,test_0069_swarm_midtime_velocity.py,test_0070_ve_stress_history_on_particles.py, a collinear-cell test intest_0067, the history-before-advection test made rank-safe. Two ranks green; nodal VEP regressions (1051, 1052) unchanged. Docs:docs/developer/subsystems/integration-point-variables.md.Underworld development team with AI support from Claude Code
🤖 Generated with Claude Code
https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G