Domain & Geometry

jno.domain holds a meshed geometry, its named regions, and the collocation points sampled on them. Build the geometry with the shape DSL (gmsh-OpenCASCADE), load a mesh file, or pass a point cloud.

import jno
from jno import shape

solid = (shape.rect(0, 0, 4, 1) - shape.disk(2, 1, 0.4)).extrude(0.6)   # a strip under a roll
d = solid.domain()                                                      # one-liner; == jno.domain(solid)

x, y, z, t = d.variable("interior")                                     # coords (+ a trailing time coord t)
xb, yb, zb, tb, nx, ny, nz = d.variable("top", normals=True, split=True)  # a named boundary + outward normals

The `shape` DSL

A shape is an immutable build-plan: make primitives, combine them with boolean operators, reshape with transforms, then hand the result to jno.domain. Every call returns a new shape.

Primitives

Each takes an optional size= (target element size near it).

Primitive	Auto-named boundaries
`shape.rect(x0, y0, x1, y1)`	`left` `right` `top` `bottom`
`shape.disk(cx, cy, r)`	`arc`
`shape.polygon(points)`	`e0 e1 … eN` (one per edge)
`shape.box(x0, y0, z0, x1, y1, z1)`	`left right top bottom front back`
`shape.cylinder(x, y, z, dx, dy, dz, r)`	`side` `top` `bottom` (any axis)
`shape.sphere(cx, cy, cz, r)`	`surface`
`Path(…).face()`	per-segment (see Contours)

interior (the area/volume) and boundary (all of it) are always present.

Combine

a - b     # cut         a | b     # fuse         a & b     # intersect

Transforms

Transform	Meaning
`.extrude(h)`	sweep a 2-D shape straight up into 3-D
`.revolve(axis_point, axis_dir, angle=2π)`	rotate a 2-D profile about an axis (x- or y-axis through the origin); a partial angle gives a wedge/half-donut
`.sweep(path)`	sweep a 2-D profile along a smooth `Path` (bent pipe)
`.array(n, step=v)` / `.array(n, about=(pt, dir), angle=2π)`	`n` fused copies — linear or polar
`.translate(v)` · `.rotate(axis_point, axis_dir, angle)`	place/orient (names are preserved)
`.fillet(radius, where=None)`	round edges; `where=f(x,y,z)` selects them by midpoint
`.sized(size)`	set the target mesh size (scalar or `f(x,y,z)`)
`.structured(n=None)`	mesh an axis-aligned rect/box as a regular lattice (see below)

extrude is a special case of sweep; revolve is a distinct rotation. Mesh size lives on the shape it describes, not on jno.domain.

Contours with arcs

Path chains line_to/arc_to segments (points are 3-D, z defaults to 0); .face() closes it into a 2-D shape. This is the general 2-D generator — a diameter plus a semicircular arc revolves into a sphere:

from jno import Path

sphere = (
    Path(0, -1)
    .line_to(0, 1, name="axis")            # the diameter, on the revolve axis
    .arc_to(0, -1, through=(1, 0), name="dome")   # a semicircular arc (3-point arc)
    .face()
    .revolve((0, 0, 0), (0, 1, 0))         # -> a sphere; the "dome" arc names its surface
)

An open Path is instead a sweep trajectory. Segments may be name='d as you draw them; the name flows onto the swept face. A sharp line_to→line_to corner is rejected (round it with arc_to).

Naming & selection

Primitives auto-name their boundaries (table above); rename, merge, or carve regions afterwards with d.tag, which accepts:

d.tag("inlet",  lambda x, y, z: x < 1e-6)                       # a coordinate predicate
d.tag("outlet", lambda x, n, name: (n[:, 0] > 0.9) & (name != "left"))  # coords + outward normal + current name

The second form selects boundary facets by position, orientation, and existing name — inclusion and exclusion in one predicate. n is the outward normal (boundary only), name each facet's current region.

A coordinate predicate works in 1-D too, where a boundary facet is a single vertex: d.tag("right_end", lambda x: x > 1 - 1e-9) names an endpoint, carries its outward normal (-1 at the left end, +1 at the right), and a u(right_end) - g term is an ordinary Dirichlet condition.

To name a region and grab its coordinates in one line, pass the predicate straight to variable (it forwards to tag, then returns the split coords):

xl, yl, zl, _ = d.variable("left", where=lambda x, y, z: x < 1e-6)   # tag AND bind in one call

Mesh density

Control element size by attaching size= to a shape — gmsh's mesh-size fields do the work (a Distance+Threshold field near a shape, a size callback for f(x,y,z)), combined by min:

shape.disk(2, 1, 0.4, size=0.02)               # fine in the band around this shape's boundary
(strip - roll).sized(0.05)                      # a global size cap for the whole shape
solid.sized(lambda x, y, z: 0.03 + 0.10 * y)    # graded: denser where the function is smaller

So "denser in a region" = give a small size= to the shape covering it, or a callable that is small there.

Cell shape — `.quad()`

Meshes are simplicial by default (triangles, tetrahedra). .quad() asks gmsh to recombine them into quadrilaterals, which works on any 2-D geometry:

d = shape.disk(0, 0, 1, size=0.1).quad().domain()      # pure quads, curved boundary and all
(plate.quad() - hole).sized(0.05)                       # survives the plan operators, like .sized()

Like size= and .curved(), this is a property of the shape rather than an argument to the solve.

3-D needs .structured(). gmsh cannot hexahedral-mesh general geometry — Recombine3DAll on a plain box returns 944 tetrahedra and no hexahedra — so hexes come from a regular lattice:

d = shape.box(0, 0, 0, 1, 1, 1).structured(n=16).quad().domain()    # 4096 hexahedra
d = shape.rect(0, 0, 1, 1).structured(n=40).quad().domain()          # 1600 quadrilaterals

.quad() and .structured() compose in either order.

Regular lattices — `.structured()`

.structured() meshes an axis-aligned rect/box as a regular grid instead of calling gmsh. Three things follow that a gmsh mesh cannot give:

hexahedra, as above — the lattice is the one 3-D plan that can be hex-meshed;
a grid descriptor on domain.grid, which is what lets jno.fdm take its assembly-free 5-/7-point stencils instead of the cotangent operator, and what a nodal field reshapes against;
exactly matched opposite faces, so a whole-domain periodic tie collapses onto one DOF rather than holding to a tolerance.

d = shape.rect(0, 0, 1, 1, size=0.1).structured().domain()   # counts from size=: 10x10 cells
d = shape.box(...).structured(n=(32, 16, 16)).domain()       # explicit, per axis

d.grid          # {"shape": (Nx, Ny[, Nz]), "spacing": (...), "origin": (...)}
u.reshape(d.grid["shape"])                                    # nodes are C-ordered

n counts cells, so a lattice has n + 1 nodes per axis and d.grid["shape"] is n + 1 — consistent with the nx/ny/nz of every other grid in jNO. A 128×128 pixel grid for a foundation model is therefore .structured(n=127). Omit n to derive it from the shape's size=.

A plan that cannot be a lattice — a CSG cut, a disk, a graded size= callable, .curved() — is refused by name rather than quietly meshed with gmsh, since the caller who then reads d.grid or expects hexes would otherwise fail somewhere else. The named faces (left/right/bottom/top, plus front/back in 3-D), boundary, interior, and any .name()/.attach() region come with the lattice as usual.

See the FEM guide's tensor-product section for what these support: volume terms and Dirichlet conditions work; surface integrals, order > 1 and h-adaptivity refuse by name.

Multi-material regions

Name each material with .name(...) and combine with + to build a multi-material domain: the sub-shapes are fragmented so element edges align exactly with every material interface (a conforming mesh), and each cell is assigned to the first region — left-to-right order is priority — whose shape contains its centroid. Each region becomes its own variable set; the outer boundary keeps its auto-names and internal interface facets are not boundary.

plate = shape.rect(0, 0, 2, 1)
core  = shape.disk(1, 0.5, 0.3)
d = (core.name("inclusion") + plate.name("matrix")).sized(0.05).domain()

d.variable("inclusion")     # the disk's cells/points (a distinct material)
d.variable("matrix")        # everything else
d.boundary_tags()           # {left, right, top, bottom, boundary} — the plate's outer edges only

+ keeps the pieces as distinct materials with a conforming interface — unlike | (fuse), which merges them into one. It composes n-ary (a + b + c), and shape.regions(inclusion=core, matrix=plate) is the equivalent keyword form. Regions may overlap: here the disk overlaps the plate, and inclusion wins inside the disk because it is listed first. Use each region tag to restrict a jno.fem term or coefficient to that material (per-region volume integration), or to sample it in a PINN. A multi-material shape is a top-level construct — call .domain() on it directly; it does not compose with boolean operators or transforms.

Region names that are not valid Python identifiers go through the dict form: shape.regions({"Quartz.1": q1, "Quartz.2": q2}). Dict order is priority order exactly as for keywords, and the two forms combine (dict entries first).

Material properties — `.attach(...)`

A region can carry its own material properties, read back off the domain as a per-region coefficient ready to drop into a weak form. d.k is exactly the per-region coefficient assembled from every region that declared a k.

This is the preferred way to express materials. It puts each property next to the region it describes, lets several properties travel together, and leaves the weak form reading in physical names (d.k, d.eps) instead of coefficients assembled elsewhere. Use attach for these when the mapping is data rather than a material — values computed elsewhere, or a case needing an explicit default= for regions that genuinely have no value.

kri = shape.polygon(v).name("Kristall").attach(k=220.0, eps=0.794)
gas = shape.polygon(w).name("Gas").attach(k=0.186, eps=1.0)
d   = (kri + gas).domain()

heat = d.k * (T.x*s.x + T.y*s.y) - d.q * s      # one equation, both materials

A value may be anything jNO treats as a coefficient: a scalar, an array, a symbolic expression, a typed view (ScalarView / VectorView / MatrixView — the view type survives, so d.eps @ u still works), or a trainable jno.np.parameter, so an attached property can be fitted or differentiated through. A plain function is also accepted and is called with the domain's spatial coordinates when the property is read (.attach(k=lambda r, z: 2.0 + 0.5*z)); it has to be deferred that way because a spatially varying coefficient is built from d.variable(...), which does not exist while the geometry plan is being written.

d.attached("eps") gives the raw {region: value} mapping instead of the coefficient, for consumers that need a dict — gap.emissivity(d.attached("eps")).

Rules worth knowing:

Apply after .name(...); like name, a later transform drops the attachment, and attaching to a shape that was never named raises.
Repeated calls merge (last wins), so properties can be built up in stages.
d.<name> raises if any region failed to declare that name, listing the ones that did not — a forgotten material surfaces at first use rather than as a region that silently conducts nothing. Use a bare d.attach(name=...) to declare the default when some regions genuinely have none. This check covers shape regions only. It reads _shape_regions, so on a mesh-file domain it never fires: a region you forgot to attach falls through to the bare-attach default rather than raising. That matters most for a coefficient whose absence changes the operator's character — a forgotten reluctivity leaves nu = 0 and with it no curl-curl term at all — so on a mesh-file domain attach every region explicitly and do not rely on this guard.
Values must agree on one view type across regions; a matrix on one region and a vector on another raises rather than silently taking whichever came last.
Properties are declared, not typed for volume regions: whether eps is used as a volume or a surface quantity is decided by the term that consumes it. Boundary tags are the exception — see d.attach(...) below, where the kind is resolved when it is declared.

A property can also be attached after the domain exists, which is the only way to attach to a domain.tag — or to a mesh-file domain, which has no shape to declare on:

d.tag("wall", lambda x, y: x < 1e-9)
d.tag("lid",  lambda x, y: y > 1 - 1e-9)
d.attach("wall", h=25.0).attach("lid", h=5.0)     # a SURFACE property, per boundary facet

Here the kind is decided once, from what the target owns on this mesh: a tag owning boundary facets is a surface quantity (d.h becomes a per-facet coefficient, see the jno.fem docs), a tag owning only cells is a volume quantity, and a tag owning both is ambiguous and raises rather than guessing — split it in two, or build the coefficient explicitly.

Interfaces between materials are auto-named by the region pair, sorted — d.variable("inclusion|matrix") gives every facet where those two materials meet (however many flat faces that spans — an E-core's air|core boundary is one tag, not one per face). Impose a coupling/flux condition there, or sample it. Interfaces are listed by d.interface_tags() and kept out of d.boundary_tags() (an interface is not the outer boundary). When an interface is topologically disjoint — two separate inclusions, or the two stacked layers of a winding — the connected pieces are additionally exposed as "a|b.0", "a|b.1", … alongside the union "a|b".

Signed distance — `.sdf(...)`

contains answers inside or not; sdf answers how far — negative inside, zero on the boundary, positive outside. It takes either an array of points or one coordinate per axis:

d = shape.sdf(points)        # (N, dim) numpy or jax array  -> an array
d = shape.sdf(x, y)          # jno Variables                -> an expression, usable in a form

The expression form is the useful one. Because the distance is exactly zero on the boundary, it turns a boundary condition that a PINN would otherwise have to learn into one that holds:

plate = jno.shape.rect(0, 0, 2, 1) - jno.shape.disk(0.7, 0.5, 0.28)
dom   = plate.domain()
x, y  = dom.variable("interior", split=True)[:2]
phi   = plate.sdf(x, y)                       # 0 on every boundary, hole included
u     = g + phi * net(x, y)                   # u == g on the boundary, for ANY network weights

No hand-written x * (1 - x) factor is needed, and none exists for a rectangle with a hole in it. The same expression gives a material that changes by region without a conforming mesh: eps = 1 + (eps_r - 1) * 0.5 * (1 - jno.np.tanh(shape.sdf(x, y) / w)).

What is exact, and what is only a bound

Distances to a single primitive — rect, box, disk, sphere, cylinder, polygon — are exact, and translating or rotating preserves them. Combining shapes (|, &, -) is exact on the boundary, which is all a boundary condition or a smoothed indicator needs, but strictly inside a combined shape the value is only a bound on the true distance. sweep and fillet have no closed form and raise; revolve is not wired yet.

Differentiable in the coordinates, not in the shape

Gradients with respect to the query point are finite everywhere, in float32 as well as float64. The shape's own radii, corners and vertices are fixed numbers, so optimising the geometry itself does not go through sdf.

Worked examples

Bolt-circle plateFilleted, drilled bracketBent pipe (sweep)Graded roll-gap

holes = shape.disk(3, 0, 0.3).array(6, about=((0, 0, 0), (0, 0, 1)))   # 6 holes in a ring
d = jno.domain((shape.rect(-5, -5, 5, 5) - holes).extrude(0.4))

part = (shape.box(0, 0, 0, 8, 4, 1) - shape.cylinder(2, 2, -1, 0, 0, 3, 0.5)) \
    .fillet(0.2, where=lambda x, y, z: z > 0.9)                          # round the top edges
d = jno.domain(part)

pipe = shape.disk(0, 0, 0.4).sweep(Path(0, 0, 0).arc_to(2, 0, 2, through=(0.6, 0, 1.4)))
d = jno.domain(pipe)

strip = shape.rect(0, 0, 4, 1, size=0.1)          # coarse in the bulk
roll  = shape.disk(2, 1, 0.4, size=0.02)          # fine near the contact arc
d = jno.domain((strip - roll).extrude(0.6))       # the carved arc is auto-named "arc"
xc, yc, zc, tc, nx, ny, nz = d.variable("arc", normals=True, split=True)   # contact points + normals

Other domain sources

Pre-meshed file — physical-group names become region tags:

d = jno.domain("part.msh")     # .msh / .vtk / .med … built anywhere (gmsh, CAD, …)

interior and boundary are derived automatically when the file does not define them, so a mesh exported from anywhere is immediately usable — no physical groups required. The boundary is topological (facets belonging to exactly one cell), which is also the only option for the many files that store no surface block at all. Names the file does define always win: a physical group called boundary keeps its own meaning, and only the missing tag is derived.

When the boundary falls into more than one connected shell, each is additionally exposed as boundary_0, boundary_1, … so an internal cavity or a second body can be addressed on its own; a single-shell part gains no numbered tags.

Two kinds of file are refused rather than half-loaded, both by name:

Mixed cell types (a mesh with both tetrahedra and hexahedra, say). jNO assembles on one element family, so it would otherwise take the first block it recognised and silently ignore the rest — measured on a real gmsh benchmark as 70 % of the domain. The error reports the blocks and how many cells would be lost.
A surface (shell) mesh — triangles in 3-D with no volume behind them, a routine CAD/STL export. Solving on a manifold is a different discretisation, not a missing setting.

Reader limits are reported as such rather than as jNO errors: an unsupported .msh version (only 2.2 and 4.1 are read) and gmsh element types meshio does not implement both say so, and say what to re-export.

Point cloud — coordinates you already have (no mesh, so no .integrate()/FD):

d = jno.domain.from_array({"interior": interior_coords, "boundary": boundary_coords})

1-D domains — an open path, exactly like every other dimension's shape: jno.Path(0, 0).line_to(1, 0).curve(size=0.01).domain() (ends named left/right). This is the form used throughout the docs; the jno.domain.line(...) shorthand still exists. jno.domain also keeps the structured grids equi_distant_rect / poseidon and the point-cloud from_array. For 2-D/3-D geometry build the shape with shape — shape.rect(...).domain(), shape.box(...).domain(), and so on. (Shapely geometries and vertex lists are also still accepted.)

Mesh-free sampling — what a PINN actually needs

shape.domain() does not mesh. A shape knows its own extent, its own membership test and its own boundary in closed form, so collocation points are drawn from the geometry directly:

d = jno.shape.box(0, 0, 0, 1, 1, 1).domain()          # no gmsh, in 1-D, 2-D or 3-D
x, y, z, t = d.variable("interior", sample=(20_000, None), split=True)

Two things differ from sampling a mesh. 20_000 means twenty thousand — there is no node set to be clipped to — and every draw is fresh, so an adaptive strategy explores the region instead of reshuffling one fixed cloud. Tagging is unchanged, and works on the boundary as well as the interior:

d.tag("hot",   lambda x, y, z: (x-.5)**2 + (y-.5)**2 + (z-.5)**2 < .04)   # a lump of interior
d.tag("inlet", lambda x, y, z: x < 1e-9)                                  # a face
bx, by, bz, t, nx, ny, nz = d.variable("inlet", sample=(500, None), normals=True, split=True)

Which of the two a predicate means is measured — jNO draws from both and sees which one it accepts — because a boundary predicate applied to an interior draw matches nothing and would look like an empty region. Boundary points land on the analytic surface, so a disk's samples lie on the circle to machine precision with exactly radial normals, not on a chord with a per-facet normal. The primitives' auto-names (left, arc, surface, …) are available before any mesh exists.

When a mesh does get built

Some things are defined on a mesh, not on the geometry: jno.fem and fem_symbols(), .integrate(), the finite-difference schemes, .points / normals_by_tag / tag_indices, and a facet predicate f(x, n, names). Reading any of them builds the mesh once, and says so:

INFO: _fem.__init__ needs a mesh; this domain was mesh-free — building it now.

Nothing you write has to change; a tag declared while mesh-free gains its mesh-derived half at that point.

Shapes with no closed form — the boundary tessellation

fillet has no analytic membership: it removes material near edges, so recursing to the child would answer for the un-filleted solid — a wrong-but-plausible mask. It is served by meshing the boundary only — the perimeter in 2-D, the surface in 3-D — which shape.tessellate() does once and caches. No volume fill runs, and contains, sample_interior and sample_boundary fall back to it automatically. The domain says which path it took:

INFO: Mesh-free domain from fillet plan (dim 3); tags [...] sample from a 372-facet boundary
tessellation.

What this costs. A facet is straight, so on a curved boundary the normal is piecewise constant and O(h) wrong — 7.5° median on a unit sphere at h=0.5, 2.2° at h=0.15, halving with h. Tighten it with .size(h). At a crease the facet normal is exact, because a facet lies wholly on one face. Membership is the polyhedron's rather than the true solid's; the two differ by the sagitta, O(h²).

sweep is not served this way. A sweep along a straight line is rewritten as an extrude and stays fully analytic, so it never needed the tessellation. A sweep along an arc is refused by name: gmsh returns that surface with its seam unsewn, and refining makes it worse — measured across h = 0.5 → 0.08 the number of edges bounding a single facet runs 7, 7, 7, 10, 13, 16, 24. Without a closed surface there is no inside to test for, so that one meshes eagerly through .build(), and the domain logs why rather than doing it silently.

Plans that stay eager

.name(...) / shape.regions(...) (region and interface tags are the mesher's conforming sub-bodies) and .structured() (already a lattice) mesh at construction as before, rather than being half-served.

translate, rotate, extrude and revolve are covered analytically — the query point is mapped into the child's frame, and a revolved boundary is drawn from the profile weighted by its swept radius (Pappus), so the inner rim of a ring is not over-sampled relative to the outer.

Sampling, time, and batching

x, y, z, t = d.variable("interior")                        # all interior mesh nodes
x, y, z, t = d.variable("interior", sample=(500, None))    # 500 sampled points
xb, yb, zb, tb, nx, ny, nz = d.variable("top", normals=True, split=True)   # boundary + outward normals

The first line — no count — is read from whether you declared a resolution. size= is that declaration: asking for a mesh of a given density and then for "the interior" unambiguously means its nodes, so it builds that mesh and hands them back, exactly as before. A convergence study that varies size therefore keeps working unchanged.

With no size= there is no declared resolution and no node set, and "no count" could equally mean collocation points or nodes — which depends on what you do afterwards. Neither is chosen for you:

ValueError: domain.variable('interior'): this domain is mesh-free and no mesh size was declared,
so the tag has no node set to hand back and no natural point count. Say which you want:
  - continuous collocation points:  variable('interior', sample=(n, None))
  - a mesh's nodes:                 give the shape a size, e.g. shape.rect(..., size=0.05),
                                    or read d.mesh first

Guessing would be silent rather than wrong-looking: finite differences over a single collocation point return a number, it is just not the one you asked for. An explicit sample=(n, None) always means mesh-free, size= or not.

variable always returns a trailing time coordinate t (a constant for steady domains), so a 3-D domain unpacks as (x, y, z, t). Time-dependent domains take time=(t0, t1, n); variable("initial") is the t=0 slice. Multiply a domain by an integer B (B * jno.domain(...), or B * shape.domain()) to replicate it across B operator-learning samples. shape.domain(**kwargs) forwards the domain arguments (time=, sample=, name=, …) — but not a constructor or mesh_size, since size lives on the shape — so a full domain is one line. d.plot("domain.png") renders the mesh, regions, and normals.

Attaching a tensor

Pass an array instead of a sampling spec to attach it as a tensor tag — the operator-learning pattern, where each batch sample carries its own field:

from jno.domain.geometries import Geometries

dom = 256 * jno.domain(constructor=Geometries.poseidon(nx=64, ny=64))
dom.variable("_f", forcing)          # (256, 64, 64, 1) — the shape you actually have

Context tensors are laid out (B, T, ...). On a steady domain the time axis is 1 by definition, so it is inserted for you: attach (B, H, W, C) and it is stored as (B, 1, H, W, C), and per-node data (B, n, k) as (B, 1, n, k). Writing the axis yourself still works and is left untouched.

This is load-bearing, not cosmetic

Without a time axis the compiler reads the first grid axis as the timestep count and hands the expression a single slice — an (4, 8, 5, 1) attach used to arrive as (5, 1), with the rest of the field discarded and no error raised.

On a time-dependent domain axis 1 is genuinely ambiguous (is it T, or a grid axis?), so a mismatch raises instead: pass T entries, or 1 to share one field across all steps.

Only tensors of rank ≥ 3 whose leading dimension is B or 1 are touched — a parameter like (B, 1, 1) or (B, F), or a shared lookup table, is left exactly as given. Multiply the domain by B before attaching, so the batch count is known.

Per-node data (B, n, k) — one value per node per sample, the operator-learning target — is stored as (B, 1, n, k) when n is the node count of one of the domain's point sets. It used to arrive as (B, k), the first node of each sample, with nothing raised. When axis 1 could be read two ways (n equals B or the timestep count) or matches no point set, the attach raises and names the layout to write instead: (B, 1, n, k) for one value per node (or any other per-sample array), (B, T, 1, k) for one value per timestep. On a time grid, (B, T, k) with no point set of T nodes is one value per sample and step, and is stored as (B, T, 1, k); a lazy source in that layout is refused (reshaping it would read it) — store it as (B, T, 1, k) yourself.

Datasets larger than memory

The same slot accepts a lazy source — anything exposing .shape and __getitem__. The handle is stored unread and sliced one batch at a time, so the dataset never has to fit:

dom.variable("_f", h5py.File("well.h5")["forcing"])     # or zarr, tensorstore, np.memmap
crux.solve(epochs=600, batchsize=32, offload_data=True)

No new dependency and no new argument — the contract is duck-typed, so jNO imports none of those libraries. offload_data=True is required, because the on-device path holds the whole array; asking for it with a lazy source raises and names the fix.

What a lazy source gives up

Indices are sorted before each gather, since h5py and zarr only accept increasing fancy indices. That reorders samples within a batch — invisible to a mean-reduced loss, but a run is not bitwise identical to one recorded before this existed.

The (B, T, ...) time axis is validated, not inserted: rewriting the layout would mean reading the source, which is the one thing it exists to avoid. Store it with the time axis, or pass an eager array.

Adaptive resampling cannot target a lazy tag — it replaces the whole point set each round. It raises rather than silently skipping.