Skip to content

peclet 1.3.0 — one partition for fluid and particles, and a distributed DEM that conserves momentum

Choose a tag to compare

@eajfpeters eajfpeters released this 27 Sep 19:38
· 120 commits to main since this release

peclet-flow 1.2.0 · peclet-dem 1.1.0 · peclet-coupling 1.1.0 ·
peclet-amr 0.2.0 · peclet-pnm 1.0.3 · peclet-voro 1.0.4 · peclet-morton 1.0.2 · peclet-geom
1.0.1 · metapackages peclet and peclet-cu13 1.3.0. peclet-halo / peclet-core 1.3.0 was
released on its own on 2026-09-26 (section below) and every member now pins its C++ headers
(PECLET_CORE_TAG v1.3.0); peclet[mpi] installs 1.3.1, which adds only the promised
DeprecationWarning on the old peclet.core.* spellings (section below; headers identical).

peclet-flow 1.2.0

Changed (named numerics changes)

  • Volumetric forces are applied where the velocity lives. A per-cell body force
    (enable_cell_force: the CFD-DEM feedback, Boussinesq buoyancy) and the drag coefficient
    (enable_drag) are now evaluated at the velocity unknown's position by one rule for every
    equation: on the staggered grid the mean of the two cells sharing the face, on the collocated grid
    the cell value. Before, the constant-density path applied the cell value half a cell off the face
    (first order), and on the variable-density, non-porous path the drag's right-hand side and its
    implicit diagonal used different values. Surface forces (pressure, viscous stress, surface
    tension) are unchanged: they are integrals over the control volume's faces, not interpolated
    cell values.
  • Variable viscosity on the staggered grid is taken on the velocity control volume's own faces.
    Each velocity component used to see the viscosity half a cell shifted. A mirror-symmetric pair of
    runs now stays symmetric to round-off (asymmetry 1e-2…3e-2 → 1e-14); in the bubble column the
    bubbles had collected at one wall.
  • Collocated variable density and surface tension. The pressure, surface-tension and body
    forces enter the implicit momentum predictor instead of a face-acceleration kick after it, and a
    balanced-force pre-projection (set_balanced_force_projection) is on by default in that
    configuration (off, with a one-time notice, when the domain has an inflow or outflow face): a
    fluid at rest under gravity with a density jump stays at rest to ~1e-12.
  • Collocated advection uses the projected, divergence-free face velocity (the
    Almgren–Bell–Colella prescription) instead of the average of the two cell velocities, so the
    advective flux is conservative. This was the last uniform-grid difference from peclet.amr.

Changed (API)

  • diagnostics.rebalance_by_weights(w) builds a weighted partition whose cuts sit on multiples of
    2^a — so the pressure multigrid can coarsen a levels in place before any redistribution —
    and returns that alignment. A code that co-decomposes with flow must build its partition from
    the same weights and this alignment (dem's migrate_to_weights(w, align=)).
  • predict_hierarchy(..., weights=) forecasts the multigrid ladder after such a rebalance and then
    returns (rows, align); without weights it returns the rows alone, exactly as before.
  • init_mpi must precede the geometry, and now raises otherwise. A geometry built first left
    each rank with a pressure hierarchy for its own block as if it were a periodic domain: every solve
    converged, but the projected velocity kept each block's divergence (max |div u| 7.5e-2 at np 2).
  • A variable viscosity with an explicitly requested velocity multigrid raises ValueError (the
    V-cycle takes a constant viscosity and used to ignore the closure silently).
  • set_superficial_velocity(enabled, axis, velocity) holds the mixture superficial velocity along
    a periodic axis.

Added

  • Coarse-level multigrid stages from core 1.2.0: the telescope's stage policy and field movement
    are core's, and a weighted level-0 partition takes a Repartition stage instead of collapsing onto
    few ranks.
  • VoF marker blocks: a curvature admissibility clip (|κ| ≤ 1/Δ_min), removal of sub-grid debris
    with its volume returned exactly to the interface, and a capillary time-step bound that
    accounts for overlapping (gas–gas) interfaces (diagnostics.set_vof_phantom_capillary_bound,
    default on). The overlap failure of colliding bubbles traced to debris curvature.
  • A faster VoF step with identical results: face-form pressure-operator storage (32 instead of 56
    bytes per cell per read), fused periodic wraps, device-resident Krylov scalars, and batched
    kernels across all local marker blocks.

Fixed

  • ε^n survives a rebalance. A redistribution whose blocks changed size reset the previous
    porosity to zero before the migrated field arrived, so the first porous projection after the move
    saw dε/dt = ε/dt (pressure off by 2e2 in coupling's rebalance test).
  • VoF markers touching a low no-slip wall no longer gain or lose volume through it (per-marker
    volume drift 3.9e-3 → 7.7e-15 over 3000 steps).
  • The pressure solve no longer discards its correction on a uniform periodic bed. When the
    right-hand side is round-off on a constant (the null space of the periodic pressure operator), the
    bottom-level CG divided 0 by 0; the NaN made the multigrid preconditioner report a non-finite
    value, and the solver zeroed the whole correction and marked the step CAPPED. The bottom CG now
    stops when p·Ap = 0, as core's GraphAMG PCG already did. Paths with p·Ap > 0 are bitwise
    unchanged; the porous fixed-bed Ergun case goes from 63 CAPPED warnings to none. Double operator
    storage had been credited with curing this; it did not.
  • Pinned to core v1.3.0.

peclet-dem 1.1.0

Changed — the distributed step

  • Momentum is conserved across rank boundaries. Each contact now has exactly one owning rank,
    and ghost contributions are summed back onto their owners at every synchronisation, so a contact
    across a block face is resolved once rather than twice or not at all. Energy created at rank faces
    is gone (a restitution test gives the single-rank kinetic energy at np 2, 4 and 8), the tangential
    (Mindlin) spring no longer flips sign when a particle migrates (momentum drift 9e-4 → 3e-8 through
    18 migrations), and the multilevel coarse cycle treats aggregates as rigid bodies, conserving
    angular momentum to ~1e-7 (was 3e-3).
  • Run-to-run reproducible at every rank count. The local particle order is canonical (by
    sending rank) instead of MPI message-arrival order, which at np ≥ 4 changed the contact
    colouring and hence the result from run to run.
  • No contact is missed across a rank face or a periodic wrap. The ghost band is the global
    contact reach, 2.1 × the largest radius anywhere (was the rank's own largest radius); with core
    1.3.0 every periodic image of a particle within the band is exchanged (13–27 missed periodic pairs
    → 0 at np 2–8).
  • The overlap projection converges to one answer. Positions are corrected by an accumulated,
    retractable push per contact: converged positions agree across rank counts to 5e-5 R (was
    1e-2 R), and clusters converge in about 3x fewer iterations.
  • The cost of all this over the previous distributed step is +2–6 % per step.

Changed — contact physics (named numerics changes)

  • Restitution follows Moreau's law: every closed contact approaching faster than the resting
    threshold targets a separation velocity of −e times its approach velocity. The previous rule gave
    a pre-separating contact the target zero, which created energy (+13 % in one step of a dense
    e = 1 cluster). The 'poisson' restitution model is unchanged.
  • Tubes and boxes are found by their circumscribed radius, so end and corner contacts are
    detected (two coaxial tubes overlapping by 0.3 used to produce no contact at all).
  • In the Hertz engine, a sphere's contact radius is its own radius, not the largest radius in the
    shape registry (in a mixture with tubes a sphere read up to 1.8x its size).
  • The overlap projection weights each contact by the translational effective mass only, matching the
    translation-only correction it applies. Spheres without spin are unchanged; non-spherical bodies
    and spinning spheres change within the solver tolerance.

Added

  • migrate_to_weights(w, align=): rebuilds exactly the aligned weighted partition flow's
    rebalance_by_weights returns, so particles and fluid own the same blocks. align=1 (default) is
    the previous weighted partition.
  • Opt-in two-way shell detection, set_shell_detection('one_way' | 'two_way') (XPBD engine): a
    thin rim or edge of one body pressed into a flat face of another, between its surface sample
    points, is seen (a tube dropped on a box falls through with 'one_way' and rests with
    'two_way'). The default 'one_way' is byte-identical to before.

Fixed

  • Two deadlocks of the distributed step at np ≥ 2–4: a rank that received no ghosts skipped the
    ghosts it owed its neighbour, and the Verlet-skin rebuild was decided per rank while the rebuild is
    collective. An empty rank no longer forces a halo rebuild on every step.
  • An explicit ghost band below the contact reach is raised to it.
  • README and PyPI page: the MPI tests run on peclet.halo; the CUDA distribution is named.

peclet-coupling 1.1.0

  • One partition for fluid and particles from construction on. CfdDem builds its first
    partition from the combined CFD + DEM cost and re-balances onto the same combined weights, with
    dem following flow's aligned partition (migrate_to_weights(w, align=)); co-location is
    asserted after every co-rebalance and at the first moving step. Before, flow and dem started on
    different partitions (48³ at np 4 put 11107 particles outside their rank's flow block). Fixed
    beds are not re-partitioned. Requires peclet-flow 1.2.0 and
    peclet-dem 1.1.0.
  • Particle radii follow their particles through every migration. Multi-rank, the driver kept
    the construction-time radius array in the pre-migration order, so after a rebalance particles
    were handed another particle's radius — wrong drag in any polydisperse moving run. The radii are
    now re-derived from dem's migrated per-particle scales.
  • Every kernel is ordered against flow's and dem's CUDA streams; the smoothing kernel iterates
    x fastest on every backend.
  • The PyPI page names the distribution it ships; pinned to core v1.3.0.

peclet-amr 0.2.0

Changed

  • Public names follow peclet.flow: Flow.set_pressure_bottom(mode=),
    Flow.set_pressure_bottom_extent(cells=), Flow.pressure_bottom_extent, predict_hierarchy
    (none shipped in 0.1.1, so no aliases), and the pressure driver is Flow.set_pressure_pcg(on).
  • Coarse-level multigrid stages are on by default. Where the distributed pressure ladder can no
    longer coarsen in place, a level merges onto sibling ranks or is repartitioned (core 1.2.0),
    falling back to a replicated tail only where nothing lifts; the coarsest level of a badly factored
    mesh is solved exactly (AMG-preconditioned CG) instead of smoothed. The multigrid now continues
    below the root brick.
  • The quadratic coarse/fine flux (cf=1) is the default on graded meshes: the solver is second
    order across a 2:1 interface. set_cf_scheme(0) reproduces pre-0.2.0 graded results (and, like
    cf=1, must precede set_solid).
  • At a 2:1 seam the convective flux reconstructs from the upwind side with level-aware samples.

Fixed

  • A cut-cell wall at a coarse octree level was 4^L too stiff (an effective viscosity μ·4^L in
    the wall cells): plane Poiseuille came out as the exact parabola shifted bodily by up to 1.4.
  • The advecting face velocity is conservative again on graded meshes (a coarse/fine correction had
    re-introduced a divergence of 7e-3).

Added

  • Flow.set_pressure_tolerance(rtol) and Flow.set_momentum_tolerance(rtol) (defaults are the
    previously hard-coded values).

peclet-pnm 1.0.3 · peclet-voro 1.0.4 · peclet-morton 1.0.2 · peclet-geom 1.0.1

No change to results or Python API.

  • voro 1.0.3 was tagged but never published: its CUDA wheel build ran the hosted runner out of
    memory (the runner itself was killed, so the fallback to fewer nvcc threads never ran). The
    CUDA wheel compiles voro's split per-subsystem sources, and the build ran six of them at once,
    each compiler peaking at up to 6.2 GB. 1.0.4 is the same code, with the CUDA wheel built from the
    split sources one compilation at a time, at one nvcc thread (identical results to the unity
    build the CPU wheels use).

  • PyPI metadata: keywords, real classifiers, project links, and a Changelog URL pointing at this
    file; one citable author identity with the ORCID.

  • READMEs (the PyPI pages) name the CUDA distribution and how to install either; voro's code sample
    calls peclet.geom, not peclet.core.geom. Every tag is gated on its landing page not naming a
    retired spelling.

  • The peclet-voro-cu13 wheel can build 3.1x faster with nvcc --threads (2268 s → 735 s on four
    cores, given the memory); the hosted release runner is capped at one thread (above). The sources
    are split into one translation unit per subsystem for development, and the wheel keeps a unity
    build so the split costs it nothing.

  • pnm, voro and geom pin core v1.3.0 (geom's vendored headers are byte-identical).

peclet-halo / peclet-core 1.3.1 — 2026-09-27 — the promised DeprecationWarning

A core-only patch released with family 1.3.0. One v* tag in peclet-core publishes
peclet-halo and the peclet-core compatibility shell at the same number. The C++ headers are
byte-identical to 1.3.0 (git diff v1.3.0 v1.3.1 -- include/ is empty), so consumers keep
PECLET_CORE_TAG v1.3.0; peclet-halo is re-released only to keep the pair at one version.

Deprecated

  • import peclet.core.geom and import peclet.core.mpi now warn — step 2 of the ladder in
    docs/CORE_BOUNDARY.md, which 1.3.0 shipped without:

    DeprecationWarning: peclet.core.geom is peclet.geom since peclet 1.2.0 — use
    from peclet import geom. peclet.core is removed in peclet 2.0.0.

    and the same for peclet.core.mpi → peclet.halo. The objects are unchanged (still the same
    objects as peclet.geom / peclet.halo); only the import warns. import peclet.core alone does
    not. The shell's pins move to peclet-geom>=1.0,<2, peclet-halo>=1.2,<2.

Tests

  • python_shell_deprecation_{geom,mpi}: python -W error::DeprecationWarning -c "import peclet.core.<old>" exits non-zero with the shell's message (an ImportError cannot pass for it),
    the default filters show the warning, and the canonical peclet.geom / peclet.halo imports are
    silent under -W error. Against the 1.3.0 shell both tests fail. Core CI now installs
    peclet-geom so the geom half runs rather than skipping.
  • The suite's own callers (flow's study gates, voro's MPI validation drivers and docs, dem's
    scene_particle docs, the quick-start notebook) use the canonical names, so no suite test sees
    the warning.

peclet-halo / peclet-core 1.3.0 — 2026-09-26 — every periodic image in the particle halo, and a faster halo build

A core-only release (the family version is unchanged). One v* tag in peclet-core publishes
peclet-halo and the peclet-core compatibility shell at the same number; the C++ headers are
versioned by the same tag, which consumers pin as PECLET_CORE_TAG. No change to the peclet.halo
Python API; the new capability is C++, for peclet.dem's distributed contact solve. Design:
dem/docs/contact_solve_framework.md §5.3 (WO-8); evidence dem/docs/contact_evidence/FOLLOWUPS.md
§4b.

Added

  • ParticleHaloTopology::build(pos, rcut, includePeriodicSelf = false, allImages = false) —
    the new, opt-in allImages sends every periodic image of a particle within rcut of another
    rank's block, one send entry per image, instead of only the image nearest that block. On a
    periodic axis the ORB leaves undecomposed (a cube at np 2 and 4, a slab split along its length)
    the nearest image is always the unshifted one, so a wrapped image inside the band was never sent
    and a pair crossing a rank face while wrapping that axis was seen by neither owner (2–3 of 5754
    pairs in dem's strong-jitter lattice at np 2 and 4). With the flag an owned index may repeat in
    one rank's send list; a receiver identifies a ghost by (owner id, shift), each shift is exactly
    0 or ±L per axis, and forward, forwardPositions and reverse (host and the Kokkos
    ParticleHalo, whose atomic reverse sums repeated entries) work per entry. It also covers a
    decomposed periodic axis whose block is narrower than twice the band.

  • FlatTopo::sendShift, parallel to sendIdx: the periodic image offset of each send entry.
    A consumer can then identify an entry by (rank, image) without re-deriving core's image
    enumeration. peclet.dem's image-aware contact ownership reads it.

Changed (faster; results unchanged)

  • ParticleMigrator::imagesWithinRcutOfBlock used to walk all 3^Dim shift combinations. It now
    prefilters the candidates per axis: a combination can qualify only if each of its per-axis gaps
    does. It returns the same set in the same order, with the same summation. This runs for every
    owned particle × every rank at every topology build. With allImages it was the dominant host
    cost of dem's distributed step: +19 % ms/step at np 8.
  • ParticleHaloTopology::build loops only over candidate ranks: those whose block comes within
    rcut of the owned particles' bounding box, images included, tested by the new
    ParticleMigrator::boxWithinRcutOfBlock. The candidates are a superset of every rank any particle
    can qualify for, so the send lists (entries and order) are unchanged. The per-particle work drops
    from O(numRanks) to O(neighbours).

Unchanged

  • With allImages off (the default) the halo is byte-identical to 1.2.0: every existing ctest,
    and a byte comparison of the full flattened topology, ghost positions, forwarded ids and reverse
    sums for two all-periodic scenes at np 1, 2, 4 and 8 against a build of the 1.2.0 headers.

Tests

  • particle_halo_images_np{1,2,4,8}: a slab periodic on every axis whose z the ORB never splits,
    against a hand-enumerated 27-image oracle — received (id, image) set exact with no duplicate,
    per-particle send-entry counts, bitwise ghost positions, reverse of image-coded values summing
    every image, the default equal to an independent reimplementation of the one-image rule, and the
    flag inert at np 1.
  • particle_halo_exchange_np{1,2,4} (Kokkos) also runs an allImages topology and requires repeated
    send entries to occur and the device forward/reverse to match the host bit for bit.
  • dem's contact and momentum batteries, 264 ctests at np 1–8, are byte-identical across the
    prefilter and the candidate-rank loop.

The peclet-core compatibility shell does not yet emit DeprecationWarnings.
docs/CORE_BOUNDARY.md ties those to the family 1.3.0; this is a core-only 1.3.0. (They
arrived in 1.3.1, above.)

peclet-halo / peclet-core 1.2.0 — 2026-09-24 — coarse-level multigrid stages, and an aligned weighted ORB

A core-only release (the family version is unchanged). Its consumers adopted it before the tag
with no core change needed: peclet-amr's pressure multigrid moves coarse levels through these
stages and turns them on by default (amr 2d3fe0e…4ee0f99), and peclet.flow's telescope
delegates to them byte-identically and gains a Repartition stage after a weighted rebalance (flow
4da7171…d907c57); both reach PyPI with those packages' next releases. One v* tag in peclet-core publishes two distributions at
the same number: peclet-halo (peclet.halo, sdist, needs MPI) and the pure-Python peclet-core
compatibility shell. The C++ headers (peclet/core/...) are versioned by the same tag, which the
consumers pin as PECLET_CORE_TAG.

No change to the peclet.halo Python API (docstrings only); the new capability is C++, for the
multigrid hierarchies of flow and peclet-amr. Design: amr/docs/amr_mg_core_boundary.md
(§4, §11); decision "Coarse-level redistribution lives in core; the hierarchies stay in the
methods" (docs/decisions/core.md).

Added

  • Coarse-level multigrid stages (include/peclet/core/decomp/): where a multigrid level goes
    when it can no longer coarsen on its own decomposition, and how its fields get there and back.
    • stage_target.hpp — chooseStageTarget(cur, levelGrid, liftable, minExtent, maxBlockCells = 0), a pure, replicated policy over the caller's lift predicate: InPlace,
      SiblingMerge (the ORB tree truncated), Repartition (a proportional ORB on np_L ranks,
      rounded up to a power of two and never above the extent cap) or Replicated. With
      maxBlockCells = 0 it is flow's telescope search verbatim (4872 ladder levels compared, all
      identical). Also shallowestLiftableMerge, repartitionTarget, largestBlockCells,
      minBlockExtent, StageKind, StageTarget, SiblingMergeChoice.
    • stage_comm.hpp — makeStageComm(parent, target) → StageComm, the stage's group and
      owners-only sub-communicators (flow's Telescope splits verbatim for a sibling merge; for a
      repartition the group is the parent itself). RAII, MPI_Finalized-guarded.
    • redistribute_topology.hpp — RedistributeTopology<Dim, T>: build the movement once from the
      two decompositions and the caller's index functors, then forward / backward every V-cycle
      with no handshake — group Gatherv/Scatterv, Allgatherv, or planned point-to-point over
      the box intersections (tags 1–10, one per level). Pure copies, bitwise.
    • gather_by_global_id.hpp — gatherByGlobalId, the id-keyed replicated gather that flow's
      GraphAMG bottom and amr's replicated tail wrote inline.
  • Aligned weighted ORB (block_decomposer.hpp): BlockDecomposer::init(numBlocks, globalSize, weights, align) balances a work field while keeping every block a multiple of align —
    coarse-first (sum onto the align-grid, weighted ORB there, refined(align) back), never snapped
    after — and is init(…, weights) bit for bit at align = 1. chooseAlignedWeighted(numBlocks, G, weights, budget = 1.05, aMax) picks the deepest align = 2^a whose weightImbalance stays
    within the budget, else today's partition.
  • A landing-pages job in release.yml that publish needs: the tag fails before upload if this
    README (the PyPI page of both distributions) names a retired spelling.

Fixed

  • redistributeGridFields no longer reads out of range when the new decomposition has fewer
    blocks than ranks (a rank beyond it only sends). Inert for every existing caller: the MPI test
    prints the same bytes at np 1/2/4/8.
  • CI on main was red from the 1.2.0 boundary split (2026-09-21, so also at the 1.1.1 tag): the
    Kokkos + Python job still imported peclet.core, which this repository no longer builds, and
    its Python ctests never ran. It imports peclet.halo now.

Documentation

  • Doxygen contracts for every new type and function: what is replicated and what is collective,
    preconditions, what throws and which throws are rank-local (and therefore fatal).
  • peclet.halo docstrings say what is collective, that MPI must already be initialised, what
    gather_ghosts does and does not return, and that reverse returns a new array. The type stub
    is regenerated to match.
  • README: an install section, the new features, measured test counts, absolute links (the
    relative ones resolved neither on PyPI nor on GitHub).

Tested

72 ctests in the plain host + MPI build, 87 with Kokkos (OpenMP), 6 Python ctests (the
state_hash byte gate passing, not skipped) — np 1–8, OMP_NUM_THREADS=2. New:
test_stage_target, test_stage_redistribute_mpi (bitwise against flow's Telescope and amr's
ReplicatedTailStage), test_stage_repartition_mpi (bitwise against redistributeGridFields),
test_aligned_weighted and _mpi.