electrostatics._mesh_plan — automatic mesh planning

tempura/electrostatics/_mesh_plan.py turns a finalized device into a mesh plan: the master lattice, the effective per-layer spacings, the outer simulation box, the realized stack bounding box, and an ordered list of refinement zones. It is internal — build_problem(...) calls it — but it is where every commensurability rule that governs the mesh lives.

The master lattice

All mesh spacings are anchored to one finest lattice, the master lattice, taken from the quantum region's resolution. A spacing \(\mathbf r = [d_x, d_y, d_z]\) is admissible only if it is an integer multiple of the master \(\mathbf m\) component-wise:

\[ \frac{\mathbf r}{\mathbf m} \in \mathbb{Z}^3 . \]

The quantum region must be the finest mesh, so any non-quantum layer requesting a finer spacing is rejected, and coarser requests are rounded up to the next master multiple (with a warning) so the lattices stay nested.

Two more rules keep the geometry commensurate with the master lattice:

• Layer thickness vs. the Z lattice. Every layer thickness must be an integer multiple of the master \(d_z\) (validated in Device.finalize()), so layer interfaces land exactly on mesh planes \(z = k\,d_z\) instead of drifting between them. A non-conforming thickness fails early, naming the layer, rather than producing an empty refinement zone deep inside build_problem().
• Device span vs. coarse spacing. Each layer's in-plane \(d_x, d_y\) must evenly tile the device length/width (_validate_span_divides_layer_resolutions), so an aligned box never clips a partial cell at the high \(+x/+y\) edge.

One globally anchored lattice

Every box used for the simulation domain, the refinement zones, and the device halo is snapped to a lattice anchored at the device corner \((-L/2, -W/2, 0)\), and every mesh pattern shares one point anchor — global_point_anchor(...), half a master cell above that corner. Because each zone spacing is an integer multiple of the master and every pattern is anchored at the same point, each coarse lattice is a structural sublattice of the master: its points are a subset of the master's. This has three consequences the old per-zone phase search had to fight for:

• there is no phase to choose — the fine grid is fixed by the master lattice alone, so it is independent of vacuum_scale and the exported plane coordinates are reproducible across domain settings;
• zones of different spacing meet without hanging nodes; and
• pescado's refine_mesh(region, pattern=...) fills each zone directly from the anchored pattern, so Tempura never enumerates candidate offsets.

The outer vacuum lattice

The simulation box extends into vacuum so the grounded boundary sits far from the device. Its spacing must stay commensurate with every resolved layer, so it is built from the per-axis least common multiple of the layer multiples,

\[ c_k = \operatorname{lcm}\!\left(\frac{r^{(1)}_k}{m_k}, \dots, \frac{r^{(L)}_k}{m_k}\right) \cdot m_k , \]

then scaled up by vacuum_resolution_scale to coarsen the far field. The vacuum extent itself is a multiple of the device size set by vacuum_scale.

Refinement zones

Zones are emitted in a deterministic order (_zone_sort_key) so meshes are reproducible: a coarse device halo around the whole stack first, then the finer quantum-region slab. Each zone records the region to refine, its target resolution, and whether it defines the protected interior that the grounded outer shell must not overlap.

Two domain modes exist. With vacuum_scale > 0 the planner builds an outer vacuum box plus a device halo (outer_vacuum mode); with vacuum_scale == 0 it wraps the stack tightly and meshes only the stack (tight_stack_bbox mode).

API

_mesh_plan

Automatic mesh planning for Tempura's active electrostatics pipeline.

global_point_anchor(device, master)

Return the single mesh-point anchor shared by every pattern.

Every automatic pattern is a Rectangular.constant lattice whose spacing is an integer multiple of master. Anchoring all of them at this one point makes each coarse lattice a structural sublattice of the master lattice: refinement zones cannot land on incommensurate phases, so the realized fine grid is independent of vacuum_scale and no hanging nodes appear where zones of different spacing meet.

The anchor is offset half a master cell from the box-aligned corner origin so that mesh points sit inside the simulation box rather than exactly on its faces, matching the master lattice's own half-cell inset.

Parameters:

Name	Type	Description	Default
device	Device	Finalized device whose XY span fixes the corner origin.	required
master	ndarray	Finest global [dx, dy, dz] lattice spacing.	required

Returns:

Type	Description
ndarray	The shared (x, y, z) point anchor for every automatic pattern.

Source code in tempura/electrostatics/_mesh_plan.py

def global_point_anchor(device: Device, master: np.ndarray) -> np.ndarray:
    """Return the single mesh-point anchor shared by every pattern.

    Every automatic pattern is a ``Rectangular.constant`` lattice whose spacing
    is an integer multiple of ``master``. Anchoring all of them at this one
    point makes each coarse lattice a structural sublattice of the master
    lattice: refinement zones cannot land on incommensurate phases, so the
    realized fine grid is independent of ``vacuum_scale`` and no hanging nodes
    appear where zones of different spacing meet.

    The anchor is offset half a *master* cell from the box-aligned corner
    origin so that mesh points sit inside the simulation box rather than exactly
    on its faces, matching the master lattice's own half-cell inset.

    Args:
        device: Finalized device whose XY span fixes the corner origin.
        master: Finest global ``[dx, dy, dz]`` lattice spacing.

    Returns:
        The shared ``(x, y, z)`` point anchor for every automatic pattern.
    """
    return _master_lattice_origin(device) + 0.5 * np.asarray(master, dtype=float)

interface_eps(master)

Return the shared-interface exclusion band, relative to the Z lattice.

Layer thicknesses are required to be integer multiples of master[2] (see :meth:Device.finalize), so every shared interface lands exactly on a mesh plane z = k * master[2]. The exclusion band only needs to be a tiny fraction of that spacing to drop points sitting on the plane while keeping the next plane up. Scaling with master keeps the band meaningful regardless of the (unspecified) physical unit of device.grid.

Parameters:

Name	Type	Description	Default
master	ndarray	Finest global [dx, dy, dz] lattice spacing.	required

Returns:

Type	Description
float	Absolute exclusion band in the same units as the mesh coordinates.

Source code in tempura/electrostatics/_mesh_plan.py

def interface_eps(master: np.ndarray) -> float:
    """Return the shared-interface exclusion band, relative to the Z lattice.

    Layer thicknesses are required to be integer multiples of ``master[2]``
    (see :meth:`Device.finalize`), so every shared interface lands exactly on a
    mesh plane ``z = k * master[2]``. The exclusion band only needs to be a tiny
    fraction of that spacing to drop points sitting on the plane while keeping
    the next plane up. Scaling with ``master`` keeps the band meaningful
    regardless of the (unspecified) physical unit of ``device.grid``.

    Args:
        master: Finest global ``[dx, dy, dz]`` lattice spacing.

    Returns:
        Absolute exclusion band in the same units as the mesh coordinates.
    """
    return 1e-6 * float(np.asarray(master, dtype=float)[2])

protected_inner_bbox(mesh_plan)

Return the bbox the outer Dirichlet shell must leave untouched.

Parameters:

Name	Type	Description	Default
mesh_plan	MeshPlan	Mapping returned by :func:resolve_mesh_plan.	required

Returns:

Type	Description
ndarray	Bounding box for the finest protected interior region when one exists,
ndarray	otherwise the full realized stack bounding box.

Source code in tempura/electrostatics/_mesh_plan.py

def protected_inner_bbox(mesh_plan: MeshPlan) -> np.ndarray:
    """Return the bbox the outer Dirichlet shell must leave untouched.

    Args:
        mesh_plan: Mapping returned by :func:`resolve_mesh_plan`.

    Returns:
        Bounding box for the finest protected interior region when one exists,
        otherwise the full realized stack bounding box.
    """
    for zone in mesh_plan["zones"]:
        if bool(zone["protects_boundary"]):
            return np.asarray(zone["region"].bbox, dtype=float).copy()
    return np.asarray(mesh_plan["stack_bbox"], dtype=float).copy()

resolve_auto_vacuum_resolution(layer_resolutions, master, *, resolution_scale=1.0)

Choose one outer-vacuum lattice commensurate with every layer.

Parameters:

Name	Type	Description	Default
layer_resolutions	Sequence[ndarray]	Effective per-layer [dx, dy, dz] mesh spacings.	required
master	ndarray	Finest global master lattice.	required

Returns:

Type	Description
ndarray	Coarse outer-vacuum spacing that remains an integer multiple of every
ndarray	resolved layer spacing.

Source code in tempura/electrostatics/_mesh_plan.py

def resolve_auto_vacuum_resolution(
    layer_resolutions: Sequence[np.ndarray],
    master: np.ndarray,
    *,
    resolution_scale: float = 1.0,
) -> np.ndarray:
    """Choose one outer-vacuum lattice commensurate with every layer.

    Args:
        layer_resolutions: Effective per-layer ``[dx, dy, dz]`` mesh spacings.
        master: Finest global master lattice.

    Returns:
        Coarse outer-vacuum spacing that remains an integer multiple of every
        resolved layer spacing.
    """
    resolution_scale = float(resolution_scale)
    if resolution_scale <= 0.0:
        raise ValueError(
            f"vacuum resolution_scale must be positive; got {resolution_scale}."
        )
    requested = resolution_scale * np.max(np.vstack(layer_resolutions), axis=0)
    common = _common_resolution_multiple(layer_resolutions, master)
    actual = np.ceil((requested / common) - 1e-12)
    actual = np.maximum(actual, 1.0) * common
    _require_integer_resolution_multiple(
        actual,
        master,
        name="vacuum_resolution",
        base_name="quantum region master lattice",
    )
    _warn_resolution_adjustment(
        name="vacuum_resolution",
        requested=requested,
        actual=actual,
        reason="rounded up to a coarse common multiple of the resolved layer spacings",
    )
    return actual

resolve_master_lattice(device, layers)

Resolve the finest global lattice from the quantum region layer resolution.

Parameters:

Name	Type	Description	Default
device	Device	Finalized device whose XY lattice defines the in-plane master spacing.	required
layers	Sequence[Layer]	Finalized device layers in stack order.	required

Returns:

Type	Description
ndarray	Finest commensurate [dx, dy, dz] lattice used as the global mesh
ndarray	anchor.

Source code in tempura/electrostatics/_mesh_plan.py

def resolve_master_lattice(device: Device, layers: Sequence[Layer]) -> np.ndarray:
    """Resolve the finest global lattice from the quantum region layer resolution.

    Args:
        device: Finalized device whose XY lattice defines the in-plane master
            spacing.
        layers: Finalized device layers in stack order.

    Returns:
        Finest commensurate ``[dx, dy, dz]`` lattice used as the global mesh
        anchor.
    """
    quantum_region_layers = [
        layer for layer in layers if isinstance(layer, QuantumRegion)
    ]
    if not quantum_region_layers:
        raise ValueError("Device must include at least one QuantumRegion layer.")

    master = np.asarray(quantum_region_layers[0].resolution, dtype=float)
    if not np.allclose(
        master[:2],
        [float(device.grid), float(device.grid)],
        atol=1e-9,
        rtol=0.0,
    ):
        raise ValueError(
            "QuantumRegion resolution must match the exact device XY grid spacing; "
            f"got {master[:2].tolist()} for device grid "
            f"{[float(device.grid), float(device.grid)]}."
        )

    for quantum_region in quantum_region_layers[1:]:
        if not np.allclose(
            master,
            np.asarray(quantum_region.resolution, dtype=float),
            atol=1e-9,
            rtol=0.0,
        ):
            raise ValueError("All QuantumRegion layers must share the same resolution.")
    return master

resolve_mesh_plan(device, layers, *, vacuum_scale, vacuum_resolution_scale=8.0, dirichlet_boundary=True)

Resolve the active automatic mesh plan for one finalized device.

Parameters:

Name	Type	Description	Default
device	Device	Finalized device whose geometry should be meshed.	required
layers	Sequence[Layer]	Finalized layers in stacking order.	required
vacuum_scale	float	Outer vacuum size multiplier relative to device size.	required
vacuum_resolution_scale	float	Coarsening factor applied to the automatic outer-vacuum lattice spacing.	8.0
dirichlet_boundary	bool	Whether build_problem() will add the optional grounded Boundary region.	True

Returns:

Type	Description
MeshPlan	Plain mapping with the resolved master lattice, effective per-layer
MeshPlan	resolutions, simulation box, realized stack bbox, and ordered
MeshPlan	automatic refinement zones.

Source code in tempura/electrostatics/_mesh_plan.py

def resolve_mesh_plan(
    device: Device,
    layers: Sequence[Layer],
    *,
    vacuum_scale: float,
    vacuum_resolution_scale: float = 8.0,
    dirichlet_boundary: bool = True,
) -> MeshPlan:
    """Resolve the active automatic mesh plan for one finalized device.

    Args:
        device: Finalized device whose geometry should be meshed.
        layers: Finalized layers in stacking order.
        vacuum_scale: Outer vacuum size multiplier relative to device size.
        vacuum_resolution_scale: Coarsening factor applied to the automatic
            outer-vacuum lattice spacing.
        dirichlet_boundary: Whether build_problem() will add the optional
            grounded Boundary region.

    Returns:
        Plain mapping with the resolved master lattice, effective per-layer
        resolutions, simulation box, realized stack bbox, and ordered
        automatic refinement zones.
    """
    vacuum_scale = float(vacuum_scale)
    if vacuum_scale < 0.0:
        raise ValueError(f"vacuum_scale must be non-negative; got {vacuum_scale}.")

    master_lattice = resolve_master_lattice(device, layers)
    layer_resolutions = {
        layer.name: _resolve_layer_resolution(layer, master_lattice) for layer in layers
    }
    _validate_span_divides_layer_resolutions(device, layer_resolutions)
    stack_bbox = _stack_bbox(layers)
    if np.isclose(vacuum_scale, 0.0, atol=0.0, rtol=0.0):
        simulation_resolution = _device_halo_resolution(
            device,
            layers,
            layer_resolutions,
        )
        simulation_box = _build_aligned_stack_box(
            device,
            stack_bbox,
            resolution=simulation_resolution,
            z_padding=0.5 * float(simulation_resolution[2])
            if dirichlet_boundary
            else 0.0,
        )
        zones = _build_mesh_zones(
            device,
            layers,
            layer_resolutions,
            simulation_box,
            stack_bbox,
            master_lattice,
            include_device_halo=False,
        )
        domain_mode = "tight_stack_bbox"
    else:
        simulation_resolution = resolve_auto_vacuum_resolution(
            list(layer_resolutions.values()),
            master_lattice,
            resolution_scale=vacuum_resolution_scale,
        )
        simulation_box = _build_aligned_vacuum_box(
            device,
            layers,
            vacuum_scale=vacuum_scale,
            resolution=simulation_resolution,
        )
        zones = _build_mesh_zones(
            device,
            layers,
            layer_resolutions,
            simulation_box,
            stack_bbox,
            master_lattice,
            include_device_halo=True,
        )
        domain_mode = "outer_vacuum"
    return {
        "domain_mode": domain_mode,
        "master_lattice": np.asarray(master_lattice, dtype=float),
        "layer_resolutions": layer_resolutions,
        "simulation_resolution": np.asarray(simulation_resolution, dtype=float),
        "simulation_box": simulation_box,
        "stack_bbox": stack_bbox,
        "zones": zones,
    }