API Reference

This page follows Tempura's importable module structure. It is a reference, not the best first read. For the full workflow, start with Triple Quantum Dot. For the shorter concept pages, see Layout, Device, Mesh, and Solve, and Theory.

from tempura import Device, Dielectric, Gate, TwoDEG, export_device_view
from tempura.electrostatics import build_problem, solve_gate_potentials
from tempura.electrostatics.solver import TFSolver, daz_to_delta
from tempura.layout import (
    crop_polygons_to_aoi,
    load_layout_polygons,
    make_aoi_bbox_from_ranges,
    polygons_to_vertices,
    rasterize_gate_vertices,
)
from tempura.layout.readers import load_layout_data
from tempura.layout.layout_pipeline import prepare_layout, build_device
from tempura.plotting import PlaneSpec, plot_problem_regions_with_mesh

Examples vs core API

The worked examples also use extract_2deg_plane(...) and extract_2deg_basis(...) from tempura.electrostatics.pescado_wrapper to reshape solved vectors for plotting. The core electrostatics entry points are still build_problem(...) and solve_gate_potentials(...).

Public spelling and implicit behavior

• Device(grid=...) is the supported public spelling for the in-plane lattice. dx and dy are still accepted as compatibility aliases.
• build_problem(...) resolves per-layer resolution=[dx, dy, dz] automatically from the finalized device. The TwoDEG sets the finest master lattice.
• build_device(...) may expand one layout source_layer into several suffixed Gate objects when that source layer contains several non-empty masks.

Module guide

• tempura: high-level layer classes plus viewer export.
• tempura.electrostatics: device construction, automatic meshing, and gate-basis solves.
• tempura.electrostatics.solver: Thomas-Fermi and lower-level linear-system helpers.
• tempura.layout: low-level polygon, AOI, and rasterization helpers.
• tempura.layout.layout_pipeline: convenience workflow from layout ROI to a finalized Device.
• tempura.plotting: mesh and region inspection helpers.
• tempura.viewer: export helpers for the browser viewer.

Package Root

Public package API for Tempura.

Device

Ordered device stack used to derive solver geometry.

Source code in tempura/electrostatics/device.py

class Device:
    """Ordered device stack used to derive solver geometry."""

    def __init__(
        self,
        *,
        length: float,
        width: float,
        grid: float | None = None,
        dx: float | None = None,
        dy: float | None = None,
    ):
        """Create a device with explicit dimensions and base XY lattice spacing."""
        self.length = _require_positive(length, name="length")
        self.width = _require_positive(width, name="width")
        self._grid = _resolve_grid_spacing(grid=grid, dx=dx, dy=dy)
        self.grid_shape = (
            _grid_count(self.width, self.grid, span_name="width", step_name="grid"),
            _grid_count(self.length, self.grid, span_name="length", step_name="grid"),
        )
        self._layers: dict[str, Layer] = {}
        self._deposition_batches: list[tuple[str, ...]] = []
        self._finalized = False

    @property
    def grid(self) -> float:
        """Return the square in-plane lattice spacing."""
        return float(self._grid)

    @property
    def dx(self) -> float:
        """Return the x lattice spacing alias for :attr:`grid`."""
        return self.grid

    @property
    def dy(self) -> float:
        """Return the y lattice spacing alias for :attr:`grid`."""
        return self.grid

    @property
    def finalized(self) -> bool:
        """Return whether :meth:`finalize` has been called successfully."""
        return self._finalized

    @property
    def layers(self) -> Mapping[str, Layer]:
        """Return the ordered layer stack as a read-only mapping."""
        return MappingProxyType(self._layers)

    def add_layer(self, layer: LayerInput) -> None:
        """Add a layer in bottom-to-top stacking order.

        Args:
            layer: Either one layer instance or one non-empty list of gates
                deposited together in the same metal batch.

        Raises:
            RuntimeError: If the device has already been finalized.
            ValueError: If names collide, a gate stencil shape is invalid, or
                two metal batches are added without an intervening non-gate
                layer.
            TypeError: If ``layer`` is not one supported layer input.
        """
        if self._finalized:
            raise RuntimeError("Cannot add layers after device.finalize().")
        batch = self._normalize_deposition_batch(layer)

        if isinstance(batch[0], Gate) and self._deposition_batches:
            previous_batch = self._deposition_batches[-1]
            previous_layer = self._layers[previous_batch[0]]
            if isinstance(previous_layer, Gate):
                raise ValueError(
                    "Gate deposition batches must be separated by a dielectric "
                    "or other non-gate layer; pass gates deposited together as "
                    "device.add_layer([...])."
                )

        names_in_batch: set[str] = set()
        resolved_resolutions: dict[int, np.ndarray] = {}
        for batch_layer in batch:
            if batch_layer.name in self._layers or batch_layer.name in names_in_batch:
                raise ValueError(
                    f"Layer name {batch_layer.name!r} is already present in the device."
                )
            names_in_batch.add(batch_layer.name)

            if batch_layer.resolution is None:
                resolved_resolutions[id(batch_layer)] = self._default_resolution_for_layer(
                    batch_layer
                )

            if isinstance(batch_layer, Gate) and batch_layer.stencil.shape != self.grid_shape:
                raise ValueError(
                    f"{batch_layer.name}.stencil shape must match device grid_shape "
                    f"{self.grid_shape}; got {batch_layer.stencil.shape}."
                )

        for batch_layer in batch:
            requested_resolution = (
                resolved_resolutions[id(batch_layer)]
                if batch_layer.resolution is None
                else batch_layer.resolution
            )
            resolution = _normalize_resolution_against_device_grid(
                requested_resolution,
                grid=self.grid,
                name=f"{batch_layer.name}.resolution",
                warn_on_z_rounding=batch_layer.resolution is not None,
            )
            object.__setattr__(
                batch_layer,
                "resolution",
                np.asarray(resolution, dtype=float),
            )
            self._layers[batch_layer.name] = batch_layer

        self._deposition_batches.append(tuple(batch_layer.name for batch_layer in batch))

    def finalize(self) -> None:
        """Realize device geometry in explicit, named stages."""
        if self._finalized:
            return

        deposition_plan = self._build_deposition_plan()
        finalization = self._realize_finalized_geometry(deposition_plan)
        self._apply_finalized_geometry(finalization)
        self._freeze_layers()
        self._finalized = True

    def _default_resolution_for_layer(self, layer: Layer) -> np.ndarray:
        """Return the default mesh resolution for one layer."""
        return _normalize_resolution_against_device_grid(
            np.array(
                [
                    float(self.grid),
                    float(self.grid),
                    min(float(self.grid), float(layer.thickness)),
                ],
                dtype=float,
            ),
            grid=self.grid,
            name=f"{layer.name}.resolution",
            warn_on_z_rounding=False,
        )

    def _normalize_deposition_batch(self, layer: LayerInput) -> list[Layer]:
        """Return one validated deposition batch from ``add_layer(...)`` input."""
        if isinstance(layer, list):
            if not layer:
                raise ValueError("Gate batches passed to add_layer() must be non-empty.")
            if not all(isinstance(batch_layer, Gate) for batch_layer in layer):
                raise TypeError(
                    "Gate batches passed to add_layer() must contain only Gate instances."
                )
            return list(layer)

        if not isinstance(layer, (Gate, Dielectric, TwoDEG)):
            raise TypeError(
                "add_layer() expects a Gate, Dielectric, TwoDEG, or non-empty list[Gate]."
            )
        return [layer]

    def _iter_deposition_batches(self) -> Iterator[list[Layer]]:
        """Yield layers grouped by deposition step."""
        for batch_names in self._deposition_batches:
            yield [self._layers[name] for name in batch_names]

    def _build_deposition_plan(self) -> tuple[DepositionBatch, ...]:
        """Return ordered deposition batches with explicit next-layer lookahead.

        Returns:
            Tuple of ``(layers, next_layer)`` pairs used by the shared
            finalization pipeline. ``next_layer`` is the next batch's leading
            layer or ``None`` for the final deposition step.
        """
        batches = list(self._iter_deposition_batches())
        return tuple(
            (
                tuple(batch),
                (batches[index + 1][0] if index + 1 < len(batches) else None),
            )
            for index, batch in enumerate(batches)
        )

    def _realize_finalized_geometry(
        self,
        deposition_plan: tuple[DepositionBatch, ...],
    ) -> DeviceFinalizationArtifacts:
        """Realize height arrays and backend shapes from one shared model.

        Args:
            deposition_plan: Ordered deposition batches with explicit
                lookahead.

        Returns:
            Mapping of layer name to finalized geometry information produced by
            one shared surface-evolution model.
        """
        surface = np.zeros(self.grid_shape, dtype=float)
        current_z0 = 0.0
        layer_geometries: DeviceFinalizationArtifacts = {}

        for batch in deposition_plan:
            layers, next_layer = batch
            lead_layer = layers[0]
            if isinstance(lead_layer, Gate):
                batch_geometries, surface, current_z0 = self._realize_gate_batch(
                    layers,
                    next_layer=next_layer,
                    surface=surface,
                    current_z0=current_z0,
                )
                layer_geometries.update(batch_geometries)
                continue

            layer_geometry, surface, current_z0 = self._realize_covering_layer(
                lead_layer,
                surface=surface,
                current_z0=current_z0,
            )
            layer_geometries[lead_layer.name] = layer_geometry

        return layer_geometries

    def _realize_gate_batch(
        self,
        gate_layers: tuple[Layer, ...],
        *,
        next_layer: Layer | None,
        surface: np.ndarray,
        current_z0: float,
    ) -> tuple[dict[str, FinalizedLayerGeometry], np.ndarray, float]:
        """Realize one metal deposition batch from the shared surface state.

        Args:
            gate_layers: Gates deposited together in one batch.
            next_layer: Leading layer in the next deposition step, when any.
            surface: Running top surface before the batch is deposited.
            current_z0: Running fallback z origin before the batch is
                deposited.

        Returns:
            Tuple of ``(geometries, next_surface, next_z0)`` for the realized
            batch.
        """
        h_before = surface.copy()
        geometries: dict[str, FinalizedLayerGeometry] = {}
        for layer in gate_layers:
            gate = layer
            mask = self._layer_mask(gate)
            z_bot = np.full_like(h_before, np.nan, dtype=float)
            z_top = np.full_like(h_before, np.nan, dtype=float)
            z_bot[mask] = h_before[mask]
            z_top[mask] = h_before[mask] + float(gate.thickness)
            geometries[gate.name] = _build_finalized_layer_geometry(
                gate,
                z_bot,
                z_top,
                device_length=self.length,
                device_width=self.width,
                fallback_z0=current_z0,
            )

        if next_layer is not None and isinstance(next_layer, Dielectric):
            next_surface = h_before.copy()
            for geometry in geometries.values():
                next_surface = np.maximum(
                    next_surface,
                    np.nan_to_num(cast(np.ndarray, geometry["z_top_arr"]), nan=-np.inf),
                )
            next_z0 = max(float(cast(float, geometry["top_z"])) for geometry in geometries.values())
        else:
            next_surface = h_before
            next_z0 = float(current_z0)

        return geometries, next_surface, next_z0

    def _realize_covering_layer(
        self,
        layer: Layer,
        *,
        surface: np.ndarray,
        current_z0: float,
    ) -> tuple[FinalizedLayerGeometry, np.ndarray, float]:
        """Realize one full-coverage layer from the shared surface state.

        Args:
            layer: Non-gate layer being deposited.
            surface: Running top surface before deposition.
            current_z0: Running fallback z origin before deposition.

        Returns:
            Tuple of ``(geometry, next_surface, next_z0)`` for the realized
            layer.
        """
        mask = self._layer_mask(layer)
        thickness = float(layer.thickness)
        base = overhang_base(surface, thickness=thickness, dx=self.grid)
        z_bot = np.full_like(surface, np.nan, dtype=float)
        z_top = np.full_like(surface, np.nan, dtype=float)
        z_bot[mask] = surface[mask]
        z_top[mask] = base[mask] + thickness
        next_surface = np.maximum(surface, np.nan_to_num(z_top, nan=-np.inf))
        geometry = _build_finalized_layer_geometry(
            layer,
            z_bot,
            z_top,
            device_length=self.length,
            device_width=self.width,
            fallback_z0=current_z0,
        )
        return geometry, next_surface, float(cast(float, geometry["top_z"]))

    def _apply_finalized_geometry(
        self,
        finalization: DeviceFinalizationArtifacts,
    ) -> None:
        """Attach realized geometry mappings to the user-facing layer specs."""
        for name, layer in self._layers.items():
            try:
                layer._set_finalized_geometry(finalization[name])
            except KeyError as exc:
                raise RuntimeError(
                    f"Missing finalized geometry for layer {name!r}."
                ) from exc

    def _freeze_layers(self) -> None:
        """Freeze all layers after the full geometry pipeline succeeds."""
        for layer in self._layers.values():
            layer._freeze()

    def _layer_mask(self, layer: Layer) -> np.ndarray:
        """Return the mask for a layer, defaulting to full coverage."""
        if isinstance(layer, Gate):
            return layer.stencil
        return np.ones(self.grid_shape, dtype=bool)

dx property

Return the x lattice spacing alias for :attr:grid.

dy property

Return the y lattice spacing alias for :attr:grid.

finalized property

Return whether :meth:finalize has been called successfully.

grid property

Return the square in-plane lattice spacing.

layers property

Return the ordered layer stack as a read-only mapping.

__init__(*, length, width, grid=None, dx=None, dy=None)

Create a device with explicit dimensions and base XY lattice spacing.

Source code in tempura/electrostatics/device.py

def __init__(
    self,
    *,
    length: float,
    width: float,
    grid: float | None = None,
    dx: float | None = None,
    dy: float | None = None,
):
    """Create a device with explicit dimensions and base XY lattice spacing."""
    self.length = _require_positive(length, name="length")
    self.width = _require_positive(width, name="width")
    self._grid = _resolve_grid_spacing(grid=grid, dx=dx, dy=dy)
    self.grid_shape = (
        _grid_count(self.width, self.grid, span_name="width", step_name="grid"),
        _grid_count(self.length, self.grid, span_name="length", step_name="grid"),
    )
    self._layers: dict[str, Layer] = {}
    self._deposition_batches: list[tuple[str, ...]] = []
    self._finalized = False

add_layer(layer)

Add a layer in bottom-to-top stacking order.

Parameters:

Name	Type	Description	Default
layer	LayerInput	Either one layer instance or one non-empty list of gates deposited together in the same metal batch.	required

Raises:

Type	Description
RuntimeError	If the device has already been finalized.
ValueError	If names collide, a gate stencil shape is invalid, or two metal batches are added without an intervening non-gate layer.
TypeError	If layer is not one supported layer input.

Source code in tempura/electrostatics/device.py

def add_layer(self, layer: LayerInput) -> None:
    """Add a layer in bottom-to-top stacking order.

    Args:
        layer: Either one layer instance or one non-empty list of gates
            deposited together in the same metal batch.

    Raises:
        RuntimeError: If the device has already been finalized.
        ValueError: If names collide, a gate stencil shape is invalid, or
            two metal batches are added without an intervening non-gate
            layer.
        TypeError: If ``layer`` is not one supported layer input.
    """
    if self._finalized:
        raise RuntimeError("Cannot add layers after device.finalize().")
    batch = self._normalize_deposition_batch(layer)

    if isinstance(batch[0], Gate) and self._deposition_batches:
        previous_batch = self._deposition_batches[-1]
        previous_layer = self._layers[previous_batch[0]]
        if isinstance(previous_layer, Gate):
            raise ValueError(
                "Gate deposition batches must be separated by a dielectric "
                "or other non-gate layer; pass gates deposited together as "
                "device.add_layer([...])."
            )

    names_in_batch: set[str] = set()
    resolved_resolutions: dict[int, np.ndarray] = {}
    for batch_layer in batch:
        if batch_layer.name in self._layers or batch_layer.name in names_in_batch:
            raise ValueError(
                f"Layer name {batch_layer.name!r} is already present in the device."
            )
        names_in_batch.add(batch_layer.name)

        if batch_layer.resolution is None:
            resolved_resolutions[id(batch_layer)] = self._default_resolution_for_layer(
                batch_layer
            )

        if isinstance(batch_layer, Gate) and batch_layer.stencil.shape != self.grid_shape:
            raise ValueError(
                f"{batch_layer.name}.stencil shape must match device grid_shape "
                f"{self.grid_shape}; got {batch_layer.stencil.shape}."
            )

    for batch_layer in batch:
        requested_resolution = (
            resolved_resolutions[id(batch_layer)]
            if batch_layer.resolution is None
            else batch_layer.resolution
        )
        resolution = _normalize_resolution_against_device_grid(
            requested_resolution,
            grid=self.grid,
            name=f"{batch_layer.name}.resolution",
            warn_on_z_rounding=batch_layer.resolution is not None,
        )
        object.__setattr__(
            batch_layer,
            "resolution",
            np.asarray(resolution, dtype=float),
        )
        self._layers[batch_layer.name] = batch_layer

    self._deposition_batches.append(tuple(batch_layer.name for batch_layer in batch))

finalize()

Realize device geometry in explicit, named stages.

Source code in tempura/electrostatics/device.py

def finalize(self) -> None:
    """Realize device geometry in explicit, named stages."""
    if self._finalized:
        return

    deposition_plan = self._build_deposition_plan()
    finalization = self._realize_finalized_geometry(deposition_plan)
    self._apply_finalized_geometry(finalization)
    self._freeze_layers()
    self._finalized = True

DeviceViewExportManifest dataclass

Paths written by :func:export_device_view.

Attributes:

Name	Type	Description
scene_name	str	Human-readable scene name stored in scene.json.
scene_json_path	Path	Path to the scene metadata file.
asset_paths	tuple[Path, ...]	Paths to binary array assets referenced by scene.json.

Source code in tempura/viewer/export.py

@dataclass(frozen=True)
class DeviceViewExportManifest:
    """Paths written by :func:`export_device_view`.

    Attributes:
        scene_name: Human-readable scene name stored in ``scene.json``.
        scene_json_path: Path to the scene metadata file.
        asset_paths: Paths to binary array assets referenced by ``scene.json``.
    """

    scene_name: str
    scene_json_path: Path
    asset_paths: tuple[Path, ...]

Dielectric dataclass

Bases: LayerBase

Dielectric layer covering the full XY device footprint.

Source code in tempura/electrostatics/device.py

@dataclass
class Dielectric(LayerBase):
    """Dielectric layer covering the full XY device footprint."""

    permittivity: float

    def __post_init__(self) -> None:
        """Validate and normalize the dielectric permittivity."""
        super().__post_init__()
        object.__setattr__(
            self,
            "permittivity",
            _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
        )

__post_init__()

Validate and normalize the dielectric permittivity.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Validate and normalize the dielectric permittivity."""
    super().__post_init__()
    object.__setattr__(
        self,
        "permittivity",
        _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
    )

Gate dataclass

Bases: LayerBase

Metallic gate layer described by a device-aligned stencil.

Source code in tempura/electrostatics/device.py

@dataclass
class Gate(LayerBase):
    """Metallic gate layer described by a device-aligned stencil."""

    stencil: np.ndarray

    def __post_init__(self) -> None:
        """Validate the gate stencil and normalize it to a boolean mask."""
        super().__post_init__()
        if self.stencil is None:
            raise ValueError("Gate requires a stencil.")
        stencil = np.asarray(self.stencil)
        if stencil.ndim != 2:
            raise ValueError(f"{self.name}.stencil must be a 2-D array; got ndim={stencil.ndim}.")
        stencil = stencil.astype(bool, copy=False)
        if not np.any(stencil):
            raise ValueError(f"{self.name}.stencil must contain at least one active cell.")
        object.__setattr__(self, "stencil", stencil)

__post_init__()

Validate the gate stencil and normalize it to a boolean mask.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Validate the gate stencil and normalize it to a boolean mask."""
    super().__post_init__()
    if self.stencil is None:
        raise ValueError("Gate requires a stencil.")
    stencil = np.asarray(self.stencil)
    if stencil.ndim != 2:
        raise ValueError(f"{self.name}.stencil must be a 2-D array; got ndim={stencil.ndim}.")
    stencil = stencil.astype(bool, copy=False)
    if not np.any(stencil):
        raise ValueError(f"{self.name}.stencil must contain at least one active cell.")
    object.__setattr__(self, "stencil", stencil)

TwoDEG dataclass

Bases: LayerBase

Two-dimensional electron gas (2DEG) layer.

Source code in tempura/electrostatics/device.py

@dataclass
class TwoDEG(LayerBase):
    """Two-dimensional electron gas (2DEG) layer."""

    permittivity: float
    boundary_condition: Literal["neumann", "helmholtz", "flexible"] = "helmholtz"
    mesh_refinement_resolution: np.ndarray | None = field(default=None, kw_only=True)

    def __post_init__(self) -> None:
        """Validate the 2DEG material parameters and boundary condition."""
        super().__post_init__()
        object.__setattr__(
            self,
            "permittivity",
            _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
        )
        valid = {"neumann", "helmholtz", "flexible"}
        if self.boundary_condition not in valid:
            raise ValueError(
                f"boundary_condition must be one of {sorted(valid)}; "
                f"got {self.boundary_condition!r}"
            )
        if self.mesh_refinement_resolution is not None:
            raise ValueError(
                f"{self.name}.mesh_refinement_resolution is no longer supported; "
                "express the finest 2DEG z mesh through resolution[2]."
            )
        object.__setattr__(self, "mesh_refinement_resolution", None)

__post_init__()

Validate the 2DEG material parameters and boundary condition.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Validate the 2DEG material parameters and boundary condition."""
    super().__post_init__()
    object.__setattr__(
        self,
        "permittivity",
        _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
    )
    valid = {"neumann", "helmholtz", "flexible"}
    if self.boundary_condition not in valid:
        raise ValueError(
            f"boundary_condition must be one of {sorted(valid)}; "
            f"got {self.boundary_condition!r}"
        )
    if self.mesh_refinement_resolution is not None:
        raise ValueError(
            f"{self.name}.mesh_refinement_resolution is no longer supported; "
            "express the finest 2DEG z mesh through resolution[2]."
        )
    object.__setattr__(self, "mesh_refinement_resolution", None)

export_device_view(device, out_dir, *, scene_name='device', colors=None, overwrite=True, verbose=False)

Export a device into files consumed by the TypeScript 3D viewer.

The exporter writes one scene.json metadata file and one or more binary array files storing layer height fields (plus gate stencils).

Parameters:

Name	Type	Description	Default
device	Device	Device to export.	required
out_dir	str | Path	Output directory for the exported scene bundle.	required
scene_name	str	Label stored in the metadata file.	'device'
colors	dict[str, str] | None	Optional per-layer color overrides keyed by layer name.	None
overwrite	bool	If False, refuse to write into a non-empty directory.	True
verbose	bool	If True, print progress messages.	False

Returns:

Type	Description
DeviceViewExportManifest	Manifest describing the files written.

Raises:

Type	Description
ValueError	If the device is empty or a color override is invalid.
FileExistsError	If overwrite is disabled and the destination exists.
RuntimeError	If the device is not finalized or finalized height arrays are missing.

Source code in tempura/viewer/export.py

def export_device_view(
    device: Device,
    out_dir: str | Path,
    *,
    scene_name: str = "device",
    colors: dict[str, str] | None = None,
    overwrite: bool = True,
    verbose: bool = False,
) -> DeviceViewExportManifest:
    """Export a device into files consumed by the TypeScript 3D viewer.

    The exporter writes one ``scene.json`` metadata file and one or more binary
    array files storing layer height fields (plus gate stencils).

    Args:
        device: Device to export.
        out_dir: Output directory for the exported scene bundle.
        scene_name: Label stored in the metadata file.
        colors: Optional per-layer color overrides keyed by layer name.
        overwrite: If ``False``, refuse to write into a non-empty directory.
        verbose: If ``True``, print progress messages.

    Returns:
        Manifest describing the files written.

    Raises:
        ValueError: If the device is empty or a color override is invalid.
        FileExistsError: If ``overwrite`` is disabled and the destination exists.
        RuntimeError: If the device is not finalized or finalized height arrays are missing.
    """

    def _log(message: str) -> None:
        """Print a progress message when verbose logging is enabled."""
        if verbose:
            print(message)

    if not device.layers:
        raise ValueError("Device must contain at least one layer before export.")

    if not device.finalized:
        raise RuntimeError("Device must be finalized before export_device_view().")

    out_path = Path(out_dir)
    if out_path.exists():
        if not out_path.is_dir():
            raise FileExistsError(f"Output path is not a directory: {out_path}")
        if not overwrite and any(out_path.iterdir()):
            raise FileExistsError(f"Output directory already exists and is not empty: {out_path}")
    else:
        out_path.mkdir(parents=True, exist_ok=True)

    arrays_dir = out_path / "arrays"
    arrays_dir.mkdir(parents=True, exist_ok=True)

    color_map = _default_layer_colors(device)
    if colors is not None:
        unknown_layers = sorted(set(colors) - set(device.layers))
        if unknown_layers:
            raise ValueError(f"Unknown layer names in colors override: {unknown_layers}")
        color_map.update(colors)

    asset_paths: list[Path] = []
    layers_payload = []
    for order, layer in enumerate(device.layers.values()):
        _log(f"[export_device_view] Exporting layer {layer.name!r}...")
        layers_payload.append(
            _export_layer(
                layer=layer,
                order=order,
                arrays_dir=arrays_dir,
                color=color_map[layer.name],
                asset_paths=asset_paths,
            )
        )

    metadata = _serialize_scene_metadata(
        device=device,
        scene_name=scene_name,
        layers_payload=layers_payload,
    )
    scene_json_path = out_path / "scene.json"
    scene_json_path.write_text(json.dumps(metadata, indent=2) + "\n")
    _log(f"[export_device_view] Wrote {scene_json_path}")

    return DeviceViewExportManifest(
        scene_name=scene_name,
        scene_json_path=scene_json_path,
        asset_paths=tuple(asset_paths),
    )

Electrostatics

Electrostatics primitives and Poisson-solver helpers.

Device

Ordered device stack used to derive solver geometry.

Source code in tempura/electrostatics/device.py

class Device:
    """Ordered device stack used to derive solver geometry."""

    def __init__(
        self,
        *,
        length: float,
        width: float,
        grid: float | None = None,
        dx: float | None = None,
        dy: float | None = None,
    ):
        """Create a device with explicit dimensions and base XY lattice spacing."""
        self.length = _require_positive(length, name="length")
        self.width = _require_positive(width, name="width")
        self._grid = _resolve_grid_spacing(grid=grid, dx=dx, dy=dy)
        self.grid_shape = (
            _grid_count(self.width, self.grid, span_name="width", step_name="grid"),
            _grid_count(self.length, self.grid, span_name="length", step_name="grid"),
        )
        self._layers: dict[str, Layer] = {}
        self._deposition_batches: list[tuple[str, ...]] = []
        self._finalized = False

    @property
    def grid(self) -> float:
        """Return the square in-plane lattice spacing."""
        return float(self._grid)

    @property
    def dx(self) -> float:
        """Return the x lattice spacing alias for :attr:`grid`."""
        return self.grid

    @property
    def dy(self) -> float:
        """Return the y lattice spacing alias for :attr:`grid`."""
        return self.grid

    @property
    def finalized(self) -> bool:
        """Return whether :meth:`finalize` has been called successfully."""
        return self._finalized

    @property
    def layers(self) -> Mapping[str, Layer]:
        """Return the ordered layer stack as a read-only mapping."""
        return MappingProxyType(self._layers)

    def add_layer(self, layer: LayerInput) -> None:
        """Add a layer in bottom-to-top stacking order.

        Args:
            layer: Either one layer instance or one non-empty list of gates
                deposited together in the same metal batch.

        Raises:
            RuntimeError: If the device has already been finalized.
            ValueError: If names collide, a gate stencil shape is invalid, or
                two metal batches are added without an intervening non-gate
                layer.
            TypeError: If ``layer`` is not one supported layer input.
        """
        if self._finalized:
            raise RuntimeError("Cannot add layers after device.finalize().")
        batch = self._normalize_deposition_batch(layer)

        if isinstance(batch[0], Gate) and self._deposition_batches:
            previous_batch = self._deposition_batches[-1]
            previous_layer = self._layers[previous_batch[0]]
            if isinstance(previous_layer, Gate):
                raise ValueError(
                    "Gate deposition batches must be separated by a dielectric "
                    "or other non-gate layer; pass gates deposited together as "
                    "device.add_layer([...])."
                )

        names_in_batch: set[str] = set()
        resolved_resolutions: dict[int, np.ndarray] = {}
        for batch_layer in batch:
            if batch_layer.name in self._layers or batch_layer.name in names_in_batch:
                raise ValueError(
                    f"Layer name {batch_layer.name!r} is already present in the device."
                )
            names_in_batch.add(batch_layer.name)

            if batch_layer.resolution is None:
                resolved_resolutions[id(batch_layer)] = self._default_resolution_for_layer(
                    batch_layer
                )

            if isinstance(batch_layer, Gate) and batch_layer.stencil.shape != self.grid_shape:
                raise ValueError(
                    f"{batch_layer.name}.stencil shape must match device grid_shape "
                    f"{self.grid_shape}; got {batch_layer.stencil.shape}."
                )

        for batch_layer in batch:
            requested_resolution = (
                resolved_resolutions[id(batch_layer)]
                if batch_layer.resolution is None
                else batch_layer.resolution
            )
            resolution = _normalize_resolution_against_device_grid(
                requested_resolution,
                grid=self.grid,
                name=f"{batch_layer.name}.resolution",
                warn_on_z_rounding=batch_layer.resolution is not None,
            )
            object.__setattr__(
                batch_layer,
                "resolution",
                np.asarray(resolution, dtype=float),
            )
            self._layers[batch_layer.name] = batch_layer

        self._deposition_batches.append(tuple(batch_layer.name for batch_layer in batch))

    def finalize(self) -> None:
        """Realize device geometry in explicit, named stages."""
        if self._finalized:
            return

        deposition_plan = self._build_deposition_plan()
        finalization = self._realize_finalized_geometry(deposition_plan)
        self._apply_finalized_geometry(finalization)
        self._freeze_layers()
        self._finalized = True

    def _default_resolution_for_layer(self, layer: Layer) -> np.ndarray:
        """Return the default mesh resolution for one layer."""
        return _normalize_resolution_against_device_grid(
            np.array(
                [
                    float(self.grid),
                    float(self.grid),
                    min(float(self.grid), float(layer.thickness)),
                ],
                dtype=float,
            ),
            grid=self.grid,
            name=f"{layer.name}.resolution",
            warn_on_z_rounding=False,
        )

    def _normalize_deposition_batch(self, layer: LayerInput) -> list[Layer]:
        """Return one validated deposition batch from ``add_layer(...)`` input."""
        if isinstance(layer, list):
            if not layer:
                raise ValueError("Gate batches passed to add_layer() must be non-empty.")
            if not all(isinstance(batch_layer, Gate) for batch_layer in layer):
                raise TypeError(
                    "Gate batches passed to add_layer() must contain only Gate instances."
                )
            return list(layer)

        if not isinstance(layer, (Gate, Dielectric, TwoDEG)):
            raise TypeError(
                "add_layer() expects a Gate, Dielectric, TwoDEG, or non-empty list[Gate]."
            )
        return [layer]

    def _iter_deposition_batches(self) -> Iterator[list[Layer]]:
        """Yield layers grouped by deposition step."""
        for batch_names in self._deposition_batches:
            yield [self._layers[name] for name in batch_names]

    def _build_deposition_plan(self) -> tuple[DepositionBatch, ...]:
        """Return ordered deposition batches with explicit next-layer lookahead.

        Returns:
            Tuple of ``(layers, next_layer)`` pairs used by the shared
            finalization pipeline. ``next_layer`` is the next batch's leading
            layer or ``None`` for the final deposition step.
        """
        batches = list(self._iter_deposition_batches())
        return tuple(
            (
                tuple(batch),
                (batches[index + 1][0] if index + 1 < len(batches) else None),
            )
            for index, batch in enumerate(batches)
        )

    def _realize_finalized_geometry(
        self,
        deposition_plan: tuple[DepositionBatch, ...],
    ) -> DeviceFinalizationArtifacts:
        """Realize height arrays and backend shapes from one shared model.

        Args:
            deposition_plan: Ordered deposition batches with explicit
                lookahead.

        Returns:
            Mapping of layer name to finalized geometry information produced by
            one shared surface-evolution model.
        """
        surface = np.zeros(self.grid_shape, dtype=float)
        current_z0 = 0.0
        layer_geometries: DeviceFinalizationArtifacts = {}

        for batch in deposition_plan:
            layers, next_layer = batch
            lead_layer = layers[0]
            if isinstance(lead_layer, Gate):
                batch_geometries, surface, current_z0 = self._realize_gate_batch(
                    layers,
                    next_layer=next_layer,
                    surface=surface,
                    current_z0=current_z0,
                )
                layer_geometries.update(batch_geometries)
                continue

            layer_geometry, surface, current_z0 = self._realize_covering_layer(
                lead_layer,
                surface=surface,
                current_z0=current_z0,
            )
            layer_geometries[lead_layer.name] = layer_geometry

        return layer_geometries

    def _realize_gate_batch(
        self,
        gate_layers: tuple[Layer, ...],
        *,
        next_layer: Layer | None,
        surface: np.ndarray,
        current_z0: float,
    ) -> tuple[dict[str, FinalizedLayerGeometry], np.ndarray, float]:
        """Realize one metal deposition batch from the shared surface state.

        Args:
            gate_layers: Gates deposited together in one batch.
            next_layer: Leading layer in the next deposition step, when any.
            surface: Running top surface before the batch is deposited.
            current_z0: Running fallback z origin before the batch is
                deposited.

        Returns:
            Tuple of ``(geometries, next_surface, next_z0)`` for the realized
            batch.
        """
        h_before = surface.copy()
        geometries: dict[str, FinalizedLayerGeometry] = {}
        for layer in gate_layers:
            gate = layer
            mask = self._layer_mask(gate)
            z_bot = np.full_like(h_before, np.nan, dtype=float)
            z_top = np.full_like(h_before, np.nan, dtype=float)
            z_bot[mask] = h_before[mask]
            z_top[mask] = h_before[mask] + float(gate.thickness)
            geometries[gate.name] = _build_finalized_layer_geometry(
                gate,
                z_bot,
                z_top,
                device_length=self.length,
                device_width=self.width,
                fallback_z0=current_z0,
            )

        if next_layer is not None and isinstance(next_layer, Dielectric):
            next_surface = h_before.copy()
            for geometry in geometries.values():
                next_surface = np.maximum(
                    next_surface,
                    np.nan_to_num(cast(np.ndarray, geometry["z_top_arr"]), nan=-np.inf),
                )
            next_z0 = max(float(cast(float, geometry["top_z"])) for geometry in geometries.values())
        else:
            next_surface = h_before
            next_z0 = float(current_z0)

        return geometries, next_surface, next_z0

    def _realize_covering_layer(
        self,
        layer: Layer,
        *,
        surface: np.ndarray,
        current_z0: float,
    ) -> tuple[FinalizedLayerGeometry, np.ndarray, float]:
        """Realize one full-coverage layer from the shared surface state.

        Args:
            layer: Non-gate layer being deposited.
            surface: Running top surface before deposition.
            current_z0: Running fallback z origin before deposition.

        Returns:
            Tuple of ``(geometry, next_surface, next_z0)`` for the realized
            layer.
        """
        mask = self._layer_mask(layer)
        thickness = float(layer.thickness)
        base = overhang_base(surface, thickness=thickness, dx=self.grid)
        z_bot = np.full_like(surface, np.nan, dtype=float)
        z_top = np.full_like(surface, np.nan, dtype=float)
        z_bot[mask] = surface[mask]
        z_top[mask] = base[mask] + thickness
        next_surface = np.maximum(surface, np.nan_to_num(z_top, nan=-np.inf))
        geometry = _build_finalized_layer_geometry(
            layer,
            z_bot,
            z_top,
            device_length=self.length,
            device_width=self.width,
            fallback_z0=current_z0,
        )
        return geometry, next_surface, float(cast(float, geometry["top_z"]))

    def _apply_finalized_geometry(
        self,
        finalization: DeviceFinalizationArtifacts,
    ) -> None:
        """Attach realized geometry mappings to the user-facing layer specs."""
        for name, layer in self._layers.items():
            try:
                layer._set_finalized_geometry(finalization[name])
            except KeyError as exc:
                raise RuntimeError(
                    f"Missing finalized geometry for layer {name!r}."
                ) from exc

    def _freeze_layers(self) -> None:
        """Freeze all layers after the full geometry pipeline succeeds."""
        for layer in self._layers.values():
            layer._freeze()

    def _layer_mask(self, layer: Layer) -> np.ndarray:
        """Return the mask for a layer, defaulting to full coverage."""
        if isinstance(layer, Gate):
            return layer.stencil
        return np.ones(self.grid_shape, dtype=bool)

dx property

Return the x lattice spacing alias for :attr:grid.

dy property

Return the y lattice spacing alias for :attr:grid.

finalized property

Return whether :meth:finalize has been called successfully.

grid property

Return the square in-plane lattice spacing.

layers property

Return the ordered layer stack as a read-only mapping.

__init__(*, length, width, grid=None, dx=None, dy=None)

Create a device with explicit dimensions and base XY lattice spacing.

Source code in tempura/electrostatics/device.py

def __init__(
    self,
    *,
    length: float,
    width: float,
    grid: float | None = None,
    dx: float | None = None,
    dy: float | None = None,
):
    """Create a device with explicit dimensions and base XY lattice spacing."""
    self.length = _require_positive(length, name="length")
    self.width = _require_positive(width, name="width")
    self._grid = _resolve_grid_spacing(grid=grid, dx=dx, dy=dy)
    self.grid_shape = (
        _grid_count(self.width, self.grid, span_name="width", step_name="grid"),
        _grid_count(self.length, self.grid, span_name="length", step_name="grid"),
    )
    self._layers: dict[str, Layer] = {}
    self._deposition_batches: list[tuple[str, ...]] = []
    self._finalized = False

add_layer(layer)

Add a layer in bottom-to-top stacking order.

Parameters:

Name	Type	Description	Default
layer	LayerInput	Either one layer instance or one non-empty list of gates deposited together in the same metal batch.	required

Raises:

Type	Description
RuntimeError	If the device has already been finalized.
ValueError	If names collide, a gate stencil shape is invalid, or two metal batches are added without an intervening non-gate layer.
TypeError	If layer is not one supported layer input.

Source code in tempura/electrostatics/device.py

def add_layer(self, layer: LayerInput) -> None:
    """Add a layer in bottom-to-top stacking order.

    Args:
        layer: Either one layer instance or one non-empty list of gates
            deposited together in the same metal batch.

    Raises:
        RuntimeError: If the device has already been finalized.
        ValueError: If names collide, a gate stencil shape is invalid, or
            two metal batches are added without an intervening non-gate
            layer.
        TypeError: If ``layer`` is not one supported layer input.
    """
    if self._finalized:
        raise RuntimeError("Cannot add layers after device.finalize().")
    batch = self._normalize_deposition_batch(layer)

    if isinstance(batch[0], Gate) and self._deposition_batches:
        previous_batch = self._deposition_batches[-1]
        previous_layer = self._layers[previous_batch[0]]
        if isinstance(previous_layer, Gate):
            raise ValueError(
                "Gate deposition batches must be separated by a dielectric "
                "or other non-gate layer; pass gates deposited together as "
                "device.add_layer([...])."
            )

    names_in_batch: set[str] = set()
    resolved_resolutions: dict[int, np.ndarray] = {}
    for batch_layer in batch:
        if batch_layer.name in self._layers or batch_layer.name in names_in_batch:
            raise ValueError(
                f"Layer name {batch_layer.name!r} is already present in the device."
            )
        names_in_batch.add(batch_layer.name)

        if batch_layer.resolution is None:
            resolved_resolutions[id(batch_layer)] = self._default_resolution_for_layer(
                batch_layer
            )

        if isinstance(batch_layer, Gate) and batch_layer.stencil.shape != self.grid_shape:
            raise ValueError(
                f"{batch_layer.name}.stencil shape must match device grid_shape "
                f"{self.grid_shape}; got {batch_layer.stencil.shape}."
            )

    for batch_layer in batch:
        requested_resolution = (
            resolved_resolutions[id(batch_layer)]
            if batch_layer.resolution is None
            else batch_layer.resolution
        )
        resolution = _normalize_resolution_against_device_grid(
            requested_resolution,
            grid=self.grid,
            name=f"{batch_layer.name}.resolution",
            warn_on_z_rounding=batch_layer.resolution is not None,
        )
        object.__setattr__(
            batch_layer,
            "resolution",
            np.asarray(resolution, dtype=float),
        )
        self._layers[batch_layer.name] = batch_layer

    self._deposition_batches.append(tuple(batch_layer.name for batch_layer in batch))

finalize()

Realize device geometry in explicit, named stages.

Source code in tempura/electrostatics/device.py

def finalize(self) -> None:
    """Realize device geometry in explicit, named stages."""
    if self._finalized:
        return

    deposition_plan = self._build_deposition_plan()
    finalization = self._realize_finalized_geometry(deposition_plan)
    self._apply_finalized_geometry(finalization)
    self._freeze_layers()
    self._finalized = True

Dielectric dataclass

Bases: LayerBase

Dielectric layer covering the full XY device footprint.

Source code in tempura/electrostatics/device.py

@dataclass
class Dielectric(LayerBase):
    """Dielectric layer covering the full XY device footprint."""

    permittivity: float

    def __post_init__(self) -> None:
        """Validate and normalize the dielectric permittivity."""
        super().__post_init__()
        object.__setattr__(
            self,
            "permittivity",
            _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
        )

__post_init__()

Validate and normalize the dielectric permittivity.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Validate and normalize the dielectric permittivity."""
    super().__post_init__()
    object.__setattr__(
        self,
        "permittivity",
        _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
    )

Gate dataclass

Bases: LayerBase

Metallic gate layer described by a device-aligned stencil.

Source code in tempura/electrostatics/device.py

@dataclass
class Gate(LayerBase):
    """Metallic gate layer described by a device-aligned stencil."""

    stencil: np.ndarray

    def __post_init__(self) -> None:
        """Validate the gate stencil and normalize it to a boolean mask."""
        super().__post_init__()
        if self.stencil is None:
            raise ValueError("Gate requires a stencil.")
        stencil = np.asarray(self.stencil)
        if stencil.ndim != 2:
            raise ValueError(f"{self.name}.stencil must be a 2-D array; got ndim={stencil.ndim}.")
        stencil = stencil.astype(bool, copy=False)
        if not np.any(stencil):
            raise ValueError(f"{self.name}.stencil must contain at least one active cell.")
        object.__setattr__(self, "stencil", stencil)

__post_init__()

Validate the gate stencil and normalize it to a boolean mask.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Validate the gate stencil and normalize it to a boolean mask."""
    super().__post_init__()
    if self.stencil is None:
        raise ValueError("Gate requires a stencil.")
    stencil = np.asarray(self.stencil)
    if stencil.ndim != 2:
        raise ValueError(f"{self.name}.stencil must be a 2-D array; got ndim={stencil.ndim}.")
    stencil = stencil.astype(bool, copy=False)
    if not np.any(stencil):
        raise ValueError(f"{self.name}.stencil must contain at least one active cell.")
    object.__setattr__(self, "stencil", stencil)

LayerBase dataclass

Shared validated inputs for device layers.

The public dataclass fields are the user-provided layer specification. Realized solver geometry is attached later by :meth:Device.finalize and exposed through read-only properties such as :attr:z_bot_arr and :attr:shape.

Source code in tempura/electrostatics/device.py

@dataclass
class LayerBase:
    """Shared validated inputs for device layers.

    The public dataclass fields are the user-provided layer specification.
    Realized solver geometry is attached later by :meth:`Device.finalize` and
    exposed through read-only properties such as :attr:`z_bot_arr` and
    :attr:`shape`.
    """

    name: str
    thickness: float
    resolution: np.ndarray | None = field(default=None, kw_only=True)

    _frozen: bool = field(init=False, default=False, repr=False)
    _finalized_geometry: FinalizedLayerGeometry | None = field(
        init=False,
        default=None,
        repr=False,
    )

    def __setattr__(self, name: str, value: object) -> None:
        """Freeze public fields once finalized while allowing private caches."""
        if getattr(self, "_frozen", False) and not name.startswith("_"):
            raise FrozenInstanceError(f"cannot assign to field {name!r}")
        object.__setattr__(self, name, value)

    def __post_init__(self) -> None:
        """Normalize the public layer specification into validated NumPy inputs."""
        object.__setattr__(self, "name", str(self.name))
        object.__setattr__(
            self,
            "thickness",
            _require_positive(self.thickness, name=f"{self.name}.thickness"),
        )
        if self.resolution is None:
            object.__setattr__(self, "resolution", None)
        else:
            object.__setattr__(
                self,
                "resolution",
                _normalize_resolution(self.resolution, name=f"{self.name}.resolution"),
            )

    @property
    def z_bot_arr(self) -> np.ndarray | None:
        """Return the finalized bottom height field, if available."""
        geometry = self._finalized_geometry
        return None if geometry is None else cast(np.ndarray, geometry["z_bot_arr"])

    @property
    def z_top_arr(self) -> np.ndarray | None:
        """Return the finalized top height field, if available."""
        geometry = self._finalized_geometry
        return None if geometry is None else cast(np.ndarray, geometry["z_top_arr"])

    @property
    def z_bot_fn(self) -> Callable[[float, float], float] | None:
        """Return the finalized bottom height sampler, if available."""
        geometry = self._finalized_geometry
        return None if geometry is None else cast(Callable[[float, float], float], geometry["z_bot_fn"])

    @property
    def z_top_fn(self) -> Callable[[float, float], float] | None:
        """Return the finalized top height sampler, if available."""
        geometry = self._finalized_geometry
        return None if geometry is None else cast(Callable[[float, float], float], geometry["z_top_fn"])

    @property
    def shape(self) -> shapes.Shape | None:
        """Return the finalized backend shape, if available."""
        geometry = self._finalized_geometry
        return None if geometry is None else cast(shapes.Shape, geometry["shape"])

    def _set_finalized_geometry(self, geometry: FinalizedLayerGeometry) -> None:
        """Attach finalized solver geometry to this layer instance."""
        object.__setattr__(self, "_finalized_geometry", geometry)

    def _freeze(self) -> None:
        """Freeze public and finalized array buffers after successful realization."""
        if isinstance(getattr(self, "stencil", None), np.ndarray):
            self.stencil.setflags(write=False)
        if isinstance(self.resolution, np.ndarray):
            self.resolution.setflags(write=False)
        mesh_refinement_resolution = getattr(self, "mesh_refinement_resolution", None)
        if isinstance(mesh_refinement_resolution, np.ndarray):
            mesh_refinement_resolution.setflags(write=False)
        geometry = self._finalized_geometry
        if geometry is not None:
            cast(np.ndarray, geometry["z_bot_arr"]).setflags(write=False)
            cast(np.ndarray, geometry["z_top_arr"]).setflags(write=False)
        object.__setattr__(self, "_frozen", True)

    def _require_finalized_geometry(self) -> FinalizedLayerGeometry:
        """Return finalized geometry or raise if finalization has not happened."""
        geometry = self._finalized_geometry
        if geometry is None:
            raise ValueError("Layer geometry is not available before device.finalize().")
        return geometry

    def _sample_height_bounds(
        self,
        x: np.ndarray,
        y: np.ndarray,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Return vectorized bottom/top height samples for world coordinates."""
        geometry = self._require_finalized_geometry()
        device_length = float(cast(float, geometry["device_length"]))
        device_width = float(cast(float, geometry["device_width"]))
        return (
            sample_height_array(
                cast(np.ndarray, geometry["z_bot_arr"]),
                device_length,
                device_width,
                x,
                y,
            ),
            sample_height_array(
                cast(np.ndarray, geometry["z_top_arr"]),
                device_length,
                device_width,
                x,
                y,
            ),
        )

shape property

Return the finalized backend shape, if available.

z_bot_arr property

Return the finalized bottom height field, if available.

z_bot_fn property

Return the finalized bottom height sampler, if available.

z_top_arr property

Return the finalized top height field, if available.

z_top_fn property

Return the finalized top height sampler, if available.

__post_init__()

Normalize the public layer specification into validated NumPy inputs.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Normalize the public layer specification into validated NumPy inputs."""
    object.__setattr__(self, "name", str(self.name))
    object.__setattr__(
        self,
        "thickness",
        _require_positive(self.thickness, name=f"{self.name}.thickness"),
    )
    if self.resolution is None:
        object.__setattr__(self, "resolution", None)
    else:
        object.__setattr__(
            self,
            "resolution",
            _normalize_resolution(self.resolution, name=f"{self.name}.resolution"),
        )

__setattr__(name, value)

Freeze public fields once finalized while allowing private caches.

Source code in tempura/electrostatics/device.py

def __setattr__(self, name: str, value: object) -> None:
    """Freeze public fields once finalized while allowing private caches."""
    if getattr(self, "_frozen", False) and not name.startswith("_"):
        raise FrozenInstanceError(f"cannot assign to field {name!r}")
    object.__setattr__(self, name, value)

TwoDEG dataclass

Bases: LayerBase

Two-dimensional electron gas (2DEG) layer.

Source code in tempura/electrostatics/device.py

@dataclass
class TwoDEG(LayerBase):
    """Two-dimensional electron gas (2DEG) layer."""

    permittivity: float
    boundary_condition: Literal["neumann", "helmholtz", "flexible"] = "helmholtz"
    mesh_refinement_resolution: np.ndarray | None = field(default=None, kw_only=True)

    def __post_init__(self) -> None:
        """Validate the 2DEG material parameters and boundary condition."""
        super().__post_init__()
        object.__setattr__(
            self,
            "permittivity",
            _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
        )
        valid = {"neumann", "helmholtz", "flexible"}
        if self.boundary_condition not in valid:
            raise ValueError(
                f"boundary_condition must be one of {sorted(valid)}; "
                f"got {self.boundary_condition!r}"
            )
        if self.mesh_refinement_resolution is not None:
            raise ValueError(
                f"{self.name}.mesh_refinement_resolution is no longer supported; "
                "express the finest 2DEG z mesh through resolution[2]."
            )
        object.__setattr__(self, "mesh_refinement_resolution", None)

__post_init__()

Validate the 2DEG material parameters and boundary condition.

Source code in tempura/electrostatics/device.py

def __post_init__(self) -> None:
    """Validate the 2DEG material parameters and boundary condition."""
    super().__post_init__()
    object.__setattr__(
        self,
        "permittivity",
        _require_positive(self.permittivity, name=f"{self.name}.permittivity"),
    )
    valid = {"neumann", "helmholtz", "flexible"}
    if self.boundary_condition not in valid:
        raise ValueError(
            f"boundary_condition must be one of {sorted(valid)}; "
            f"got {self.boundary_condition!r}"
        )
    if self.mesh_refinement_resolution is not None:
        raise ValueError(
            f"{self.name}.mesh_refinement_resolution is no longer supported; "
            "express the finest 2DEG z mesh through resolution[2]."
        )
    object.__setattr__(self, "mesh_refinement_resolution", None)

build_problem(device, vacuum_scale=8.0, *, verbose=False)

Build a Poisson problem for a layered device.

Parameters:

Name	Type	Description	Default
device	Device	Finalized device whose electrostatic stack should be meshed.	required
vacuum_scale	float	Outer vacuum size multiplier relative to the device span.	8.0
verbose	bool	Whether to emit progress logging during assembly.	False

Returns:

Type	Description
ProblemBuilder	Tuple of (problem_builder, region_shapes, gate_names) compatible
RegionMap	with Tempura's historical public API.

Source code in tempura/electrostatics/pescado_wrapper.py

def build_problem(
    device: Device,
    vacuum_scale: float = 8.0,
    *,
    verbose: bool = False,
) -> tuple[ProblemBuilder, RegionMap, list[str]]:
    """Build a Poisson problem for a layered device.

    Args:
        device: Finalized device whose electrostatic stack should be meshed.
        vacuum_scale: Outer vacuum size multiplier relative to the device span.
        verbose: Whether to emit progress logging during assembly.

    Returns:
        Tuple of ``(problem_builder, region_shapes, gate_names)`` compatible
        with Tempura's historical public API.
    """
    problem_builder, region_shapes, gate_names, _ = build_problem_components(
        device,
        vacuum_scale=vacuum_scale,
        verbose=verbose,
    )
    return problem_builder, region_shapes, gate_names

build_rhs_vectors(region_inds, gate_names, *, size=None, charge=None, dtype=DTYPE)

Build the multi-RHS matrix for the ordered gate basis solves.

Source code in tempura/electrostatics/solver.py

def build_rhs_vectors(
    region_inds: IndexMap,
    gate_names: list[str],
    *,
    size: int | None = None,
    charge: np.ndarray | None = None,
    dtype: DTypeLike = DTYPE,
) -> sp.csc_matrix | np.ndarray:
    """Build the multi-RHS matrix for the ordered gate basis solves."""
    gate_inds = _ordered_gate_inds(region_inds, gate_names)
    n_gates = len(gate_inds)
    n_rows = _region_system_size(region_inds, gate_inds) if size is None else size
    if n_gates == 0:
        return sp.csc_matrix((n_rows, 0), dtype=np.dtype(dtype))

    gate_sizes = np.fromiter(
        (inds.size for inds in gate_inds), dtype=int, count=n_gates
    )
    gate_row_inds = (
        np.concatenate(gate_inds) if gate_sizes.sum() else np.empty(0, dtype=int)
    )
    gate_col_inds = np.repeat(np.arange(n_gates), gate_sizes)

    # Without distributed charge, each RHS column is just a sparse gate basis
    # vector, so keep the representation sparse as long as possible.
    if charge is None:
        return sp.csc_matrix(
            (
                np.ones(gate_row_inds.size, dtype=np.dtype(dtype)),
                (gate_row_inds, gate_col_inds),
            ),
            shape=(n_rows, n_gates),
        )

    charge_matrix = np.asarray(charge, dtype=np.dtype(dtype))
    if charge_matrix.ndim == 1:
        charge_matrix = charge_matrix[:, None]
    if charge_matrix.shape[0] != n_rows:
        raise ValueError(
            f"charge must have {n_rows} rows; got {charge_matrix.shape[0]}"
        )
    if charge_matrix.shape[1] == 1 and n_gates != 1:
        charge_matrix = np.repeat(charge_matrix, n_gates, axis=1)
    elif charge_matrix.shape[1] != n_gates:
        raise ValueError(
            f"charge must have {n_gates} columns; got {charge_matrix.shape[1]}"
        )

    # Charge terms populate interior rows for every solve, so the dense path is
    # simpler than building a mixed sparse/dense block structure.
    rhs_matrix = np.zeros((n_rows, n_gates), dtype=np.dtype(dtype), order="F")
    interior_inds = np.asarray(region_inds["Interior"], dtype=int)
    if interior_inds.size:
        rhs_matrix[interior_inds, :] = charge_matrix[interior_inds, :]
    if gate_row_inds.size:
        rhs_matrix[gate_row_inds, gate_col_inds] = 1.0
    return rhs_matrix

get_region_inds(pescado_problem_finalized, gate_shapes, boundary_shape, *, boundary_name='Boundary')

Return indices for interior, gates, and the outer boundary.

Source code in tempura/electrostatics/_pescado_solve.py

def get_region_inds(
    pescado_problem_finalized: Any,
    gate_shapes: RegionMap,
    boundary_shape: shapes.Shape,
    *,
    boundary_name: str = "Boundary",
) -> IndexMap:
    """Return indices for interior, gates, and the outer boundary."""
    region_inds: IndexMap = {}
    region_inds["Interior"] = np.concatenate(
        (
            np.asarray(pescado_problem_finalized.neumann_indices, dtype=int),
            np.asarray(pescado_problem_finalized.helmholtz_indices, dtype=int),
        )
    )
    for gate_name, gate_shape in gate_shapes.items():
        region_inds[gate_name] = pescado_problem_finalized.points_inside(gate_shape)

    region_inds[boundary_name] = pescado_problem_finalized.points_inside(boundary_shape)
    return region_inds

setup_solver(matrix, *, dtype=DTYPE, blr=BLR, eps_blr=EPS_BLR, symmetric=False)

Create and factor a MUMPS solver for a sparse linear system.

Parameters:

Name	Type	Description	Default
matrix	Any	Real-valued system matrix to factor.	required
dtype	DTypeLike	Real numeric dtype used to cast matrix before factorization. Tempura's electrostatics solve path does not preserve imaginary components.	DTYPE
blr	bool	Whether to enable MUMPS block low-rank factorization.	BLR
eps_blr	float	BLR compression threshold used when blr is enabled.	EPS_BLR
symmetric	bool	Whether matrix should be factorized as symmetric.	False

Returns:

Type	Description
Any	Factored MUMPS context that can be reused across solves.

Raises:

Type	Description
RuntimeError	If python-mumps is not available.

Source code in tempura/electrostatics/solver.py

def setup_solver(
    matrix: Any,
    *,
    dtype: DTypeLike = DTYPE,
    blr: bool = BLR,
    eps_blr: float = EPS_BLR,
    symmetric: bool = False,
) -> Any:
    """Create and factor a MUMPS solver for a sparse linear system.

    Args:
        matrix: Real-valued system matrix to factor.
        dtype: Real numeric dtype used to cast ``matrix`` before factorization.
            Tempura's electrostatics solve path does not preserve imaginary
            components.
        blr: Whether to enable MUMPS block low-rank factorization.
        eps_blr: BLR compression threshold used when ``blr`` is enabled.
        symmetric: Whether ``matrix`` should be factorized as symmetric.

    Returns:
        Factored MUMPS context that can be reused across solves.

    Raises:
        RuntimeError: If ``python-mumps`` is not available.
    """
    solver = _require_mumps().Context()
    _factorize_solver(
        solver,
        matrix,
        dtype=dtype,
        blr=blr,
        eps_blr=eps_blr,
        symmetric=symmetric,
    )
    return solver

solve_gate_potentials(problem_builder, region_shapes, gate_names, *, charge=None, dtype=DTYPE, blr=BLR, eps_blr=EPS_BLR, rhs_block_size=8, save_2deg_potential=None, twodeg_plane_selection='midpoint', verbose=False)

Solve the gate basis problem with one reused factorization.

Parameters:

Name	Type	Description	Default
problem_builder	ProblemBuilder	Initialized Pescado problem builder returned by :func:build_problem.	required
region_shapes	RegionMap	Mapping of region names to realized region shapes.	required
gate_names	list[str]	Ordered gate names to solve.	required
charge	ndarray | None	Optional distributed charge columns to include in the RHS.	None
dtype	DTypeLike	Numeric dtype used during factorization and solve.	DTYPE
blr	bool	Whether to enable MUMPS block low-rank factorization.	BLR
eps_blr	float	BLR compression threshold when blr is enabled.	EPS_BLR
rhs_block_size	int	Number of gate basis columns to solve per linear solve.	8
save_2deg_potential	str | Path | None	Optional output directory for static 2DEG export.	None
twodeg_plane_selection	PlaneSelection	Which realized 2DEG z plane to export when the region spans multiple z coordinates.	'midpoint'
verbose	bool	Whether to emit progress logging during solving.	False

Returns:

Type	Description
PotentialMap	Tuple of (basis_potentials, finalized_problem) compatible with the
Any	historical public API.

Source code in tempura/electrostatics/pescado_wrapper.py

def solve_gate_potentials(
    problem_builder: ProblemBuilder,
    region_shapes: RegionMap,
    gate_names: list[str],
    *,
    charge: np.ndarray | None = None,
    dtype: DTypeLike = DTYPE,
    blr: bool = BLR,
    eps_blr: float = EPS_BLR,
    rhs_block_size: int = 8,
    save_2deg_potential: str | Path | None = None,
    twodeg_plane_selection: PlaneSelection = "midpoint",
    verbose: bool = False,
) -> tuple[PotentialMap, Any]:
    """Solve the gate basis problem with one reused factorization.

    Args:
        problem_builder: Initialized Pescado problem builder returned by
            :func:`build_problem`.
        region_shapes: Mapping of region names to realized region shapes.
        gate_names: Ordered gate names to solve.
        charge: Optional distributed charge columns to include in the RHS.
        dtype: Numeric dtype used during factorization and solve.
        blr: Whether to enable MUMPS block low-rank factorization.
        eps_blr: BLR compression threshold when ``blr`` is enabled.
        rhs_block_size: Number of gate basis columns to solve per linear solve.
        save_2deg_potential: Optional output directory for static 2DEG export.
        twodeg_plane_selection: Which realized 2DEG z plane to export when the
            region spans multiple z coordinates.
        verbose: Whether to emit progress logging during solving.

    Returns:
        Tuple of ``(basis_potentials, finalized_problem)`` compatible with the
        historical public API.
    """
    basis_potentials, problem, _ = solve_gate_potentials_components(
        problem_builder,
        region_shapes,
        gate_names,
        charge=charge,
        dtype=dtype,
        blr=blr,
        eps_blr=eps_blr,
        rhs_block_size=rhs_block_size,
        save_2deg_potential=save_2deg_potential,
        twodeg_plane_selection=twodeg_plane_selection,
        verbose=verbose,
    )
    return basis_potentials, problem

solve_linear_system(matrix, rhs, *, solver=None, dtype=DTYPE, blr=BLR, eps_blr=EPS_BLR, symmetric=False)

Solve matrix x = rhs with a factored MUMPS solver.

Parameters:

Name	Type	Description	Default
matrix	Any	Real-valued system matrix to solve.	required
rhs	Any	Dense or sparse right-hand side array.	required
solver	Any | None	Optional existing factored MUMPS context to reuse.	None
dtype	DTypeLike	Real numeric dtype used to cast the matrix and RHS. Tempura's electrostatics solve path does not preserve imaginary components.	DTYPE
blr	bool	Whether to enable MUMPS block low-rank factorization.	BLR
eps_blr	float	BLR compression threshold used when blr is enabled.	EPS_BLR
symmetric	bool	Whether matrix should be solved as symmetric.	False

Returns:

Type	Description
tuple[ndarray, Any]	Tuple of (solution, solver).

Raises:

Type	Description
RuntimeError	If python-mumps is not available and solver is not provided.

Source code in tempura/electrostatics/solver.py

def solve_linear_system(
    matrix: Any,
    rhs: Any,
    *,
    solver: Any | None = None,
    dtype: DTypeLike = DTYPE,
    blr: bool = BLR,
    eps_blr: float = EPS_BLR,
    symmetric: bool = False,
) -> tuple[np.ndarray, Any]:
    """Solve ``matrix x = rhs`` with a factored MUMPS solver.

    Args:
        matrix: Real-valued system matrix to solve.
        rhs: Dense or sparse right-hand side array.
        solver: Optional existing factored MUMPS context to reuse.
        dtype: Real numeric dtype used to cast the matrix and RHS. Tempura's
            electrostatics solve path does not preserve imaginary components.
        blr: Whether to enable MUMPS block low-rank factorization.
        eps_blr: BLR compression threshold used when ``blr`` is enabled.
        symmetric: Whether ``matrix`` should be solved as symmetric.

    Returns:
        Tuple of ``(solution, solver)``.

    Raises:
        RuntimeError: If ``python-mumps`` is not available and ``solver`` is
            not provided.
    """
    if solver is None:
        solver = _require_mumps().Context()

    _factorize_solver(
        solver,
        matrix,
        dtype=dtype,
        blr=blr,
        eps_blr=eps_blr,
        symmetric=symmetric,
    )

    if sp.issparse(rhs):
        rhs_fortran: Any = sp.csc_matrix(rhs, dtype=np.dtype(dtype))
    else:
        rhs_fortran = np.asfortranarray(np.asarray(rhs, dtype=np.dtype(dtype)))
    solution = np.asarray(solver.solve(rhs_fortran))
    if np.issubdtype(np.dtype(dtype), np.complexfloating):
        solution = solution.real
    return solution, solver

Return indices for interior, gates, and the outer boundary.

Source code in tempura/electrostatics/_pescado_solve.py

def get_region_inds(
    pescado_problem_finalized: Any,
    gate_shapes: RegionMap,
    boundary_shape: shapes.Shape,
    *,
    boundary_name: str = "Boundary",
) -> IndexMap:
    """Return indices for interior, gates, and the outer boundary."""
    region_inds: IndexMap = {}
    region_inds["Interior"] = np.concatenate(
        (
            np.asarray(pescado_problem_finalized.neumann_indices, dtype=int),
            np.asarray(pescado_problem_finalized.helmholtz_indices, dtype=int),
        )
    )
    for gate_name, gate_shape in gate_shapes.items():
        region_inds[gate_name] = pescado_problem_finalized.points_inside(gate_shape)

    region_inds[boundary_name] = pescado_problem_finalized.points_inside(boundary_shape)
    return region_inds

Return the rectangular 2DEG plane extracted from problem.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem exposing coordinates and points_inside.	required
region_shapes	Mapping[str, Any]	Mapping of region names to shapes. Must contain "2deg".	required
plane_selection	PlaneSelection	Which z plane to export when the realized 2DEG spans multiple z coordinates. Use "lower", "midpoint", "upper", or an explicit z value.	'midpoint'

Returns:

Type	Description
TwoDEGPlaneData	Sorted rectangular plane data for the 2DEG region.

Raises:

Type	Description
ValueError	If the 2DEG region is missing, empty, or does not form a rectangular plane.

Source code in tempura/electrostatics/_twodeg.py

def extract_2deg_plane(
    problem: Any,
    region_shapes: Mapping[str, Any],
    *,
    plane_selection: PlaneSelection = "midpoint",
) -> TwoDEGPlaneData:
    """Return the rectangular 2DEG plane extracted from ``problem``.

    Args:
        problem: Finalized Pescado problem exposing ``coordinates`` and
            ``points_inside``.
        region_shapes: Mapping of region names to shapes. Must contain ``"2deg"``.
        plane_selection: Which z plane to export when the realized 2DEG spans
            multiple z coordinates. Use ``"lower"``, ``"midpoint"``,
            ``"upper"``, or an explicit z value.

    Returns:
        Sorted rectangular plane data for the 2DEG region.

    Raises:
        ValueError: If the 2DEG region is missing, empty, or does not form a
            rectangular plane.
    """
    try:
        region_shape = region_shapes["2deg"]
    except KeyError as exc:
        raise ValueError("2DEG export requires a region_shapes['2deg'] entry.") from exc

    indices_2deg = np.asarray(problem.points_inside(region_shape), dtype=int)
    coords_2deg = np.asarray(problem.coordinates[indices_2deg], dtype=float)
    if coords_2deg.size == 0:
        raise ValueError("No solver coordinates were found inside the 2DEG region.")

    z_vals = np.unique(np.round(coords_2deg[:, 2], 12))
    if len(z_vals) > 1:
        z_target = _resolve_plane_z_value(z_vals, plane_selection=plane_selection)
        plane_mask = np.isclose(coords_2deg[:, 2], z_target, atol=1e-12, rtol=0.0)
        if not np.any(plane_mask):
            raise ValueError(
                "The requested 2DEG plane does not match any realized z coordinates; "
                f"requested={z_target}, available={z_vals.tolist()}."
            )
        coords_2deg = coords_2deg[plane_mask]
        indices_2deg = indices_2deg[plane_mask]

    order = np.lexsort((coords_2deg[:, 0], coords_2deg[:, 1]))
    coords_2deg = coords_2deg[order]
    sorted_indices_2deg = indices_2deg[order]

    x_values = np.unique(np.round(coords_2deg[:, 0], 12))
    y_values = np.unique(np.round(coords_2deg[:, 1], 12))
    nx = len(x_values)
    ny = len(y_values)
    if coords_2deg.shape[0] != nx * ny:
        raise ValueError("The extracted 2DEG coordinates do not form a rectangular grid.")

    x_grid = coords_2deg[:, 0].reshape(ny, nx)
    y_grid = coords_2deg[:, 1].reshape(ny, nx)
    if not np.allclose(x_grid, np.broadcast_to(x_values, (ny, nx))):
        raise ValueError("The extracted 2DEG x coordinates do not form a rectangular grid.")
    if not np.allclose(y_grid, np.broadcast_to(y_values[:, None], (ny, nx))):
        raise ValueError("The extracted 2DEG y coordinates do not form a rectangular grid.")

    return TwoDEGPlaneData(
        sorted_indices=sorted_indices_2deg,
        coords=coords_2deg,
        x_values=np.asarray(x_values, dtype=float),
        y_values=np.asarray(y_values, dtype=float),
        z_value=float(coords_2deg[0, 2]),
        nx=nx,
        ny=ny,
    )

Reshape solved gate basis vectors onto a rectangular 2DEG plane.

Parameters:

Name	Type	Description	Default
basis_potentials	Mapping[str, ndarray]	Mapping of gate name to solved potential vector.	required
gate_names	list[str]	Gate order to preserve in the returned mapping.	required
plane	TwoDEGPlaneData	Extracted 2DEG plane metadata and sorted solver indices.	required

Returns:

Type	Description
dict[str, ndarray]	Mapping of gate name to (ny, nx) 2DEG basis potential arrays.

Source code in tempura/electrostatics/_twodeg.py

def extract_2deg_basis(
    basis_potentials: Mapping[str, np.ndarray],
    gate_names: list[str],
    plane: TwoDEGPlaneData,
) -> dict[str, np.ndarray]:
    """Reshape solved gate basis vectors onto a rectangular 2DEG plane.

    Args:
        basis_potentials: Mapping of gate name to solved potential vector.
        gate_names: Gate order to preserve in the returned mapping.
        plane: Extracted 2DEG plane metadata and sorted solver indices.

    Returns:
        Mapping of gate name to ``(ny, nx)`` 2DEG basis potential arrays.
    """
    return {
        gate_name: np.asarray(basis_potentials[gate_name])[plane.sorted_indices].reshape(
            plane.ny,
            plane.nx,
        )
        for gate_name in gate_names
    }

Thomas-Fermi wrapper for self-consistent 2DEG solves.

Source code in tempura/electrostatics/solver.py

class TFSolver:
    """Thomas-Fermi wrapper for self-consistent 2DEG solves."""

    def __init__(
        self,
        pd: Any,
        twodeg_vol: np.ndarray,
        twodeg_height: float,
        twodeg_idx: np.ndarray | None = None,
        daz: float | None = None,
        me_eff: float = 0.022,
        *,
        verbose: bool = False,
    ) -> None:
        self.pd = pd
        self.verbose = bool(verbose)
        self.twodeg_idx = (
            self.pd.points(name="2deg") if twodeg_idx is None else twodeg_idx
        )
        self.twodeg_idx = np.asarray(self.twodeg_idx, dtype=int)
        if daz is not None:
            delta = daz_to_delta(daz, me_eff=me_eff)
        else:
            delta = 0.05
        self.delta = float(delta)
        if self.verbose:
            print(f"Using delta = {self.delta:.3f} eV for the ILDOS calculation.")
        ildos_arr_list, sites_ildos = twodeg_ILDOS(
            pd=self.pd,
            thickness=twodeg_height,
            twodeg_vol=twodeg_vol,
            delta=self.delta,
            me_eff=me_eff,
            twodeg_idx=self.twodeg_idx,
        )
        tf_plinear = tf_backend.thomas_fermi(
            ildos=ildos_arr_list,
            sites_ildos=sites_ildos,
            poisson_problem=self.pd,
            e_f=0.0,
        )
        self.solver = tf_plinear

    def solve(
        self,
        voltages: dict[str, float],
        initial_guess: SparseVector | None = None,
        max_ite: int = 15,
    ) -> None:
        """Solve the self-consistent problem for the provided region voltages."""
        sv_gate = self.pd.sparse_vector(val=0, name="Boundary")
        for key, value in voltages.items():
            sv_gate.extend(self.pd.sparse_vector(val=float(value), name=key))
        sv_charge = self.pd.sparse_vector(val=0, name="Boundary")
        poisson_input = {"voltage": sv_gate, "charge": sv_charge}
        if initial_guess is None:
            initial_guess = SparseVector(
                values=np.zeros(len(self.twodeg_idx)),
                indices=self.twodeg_idx,
                dtype=float,
            )

        self.solver.solve(
            poisson_input=poisson_input, initial_guess=initial_guess, max_ite=max_ite
        )

    def get_potential(self) -> np.ndarray:
        """Return the latest self-consistent chemical potential field."""
        return self.solver.chemical_potential(iteration=-1)

    def get_charge(self) -> np.ndarray:
        """Return the latest self-consistent quantum charge field."""
        return self.solver.quantum_charge(iteration=-1)

get_charge()

Return the latest self-consistent quantum charge field.

Source code in tempura/electrostatics/solver.py

def get_charge(self) -> np.ndarray:
    """Return the latest self-consistent quantum charge field."""
    return self.solver.quantum_charge(iteration=-1)

get_potential()

Return the latest self-consistent chemical potential field.

Source code in tempura/electrostatics/solver.py

def get_potential(self) -> np.ndarray:
    """Return the latest self-consistent chemical potential field."""
    return self.solver.chemical_potential(iteration=-1)

solve(voltages, initial_guess=None, max_ite=15)

Solve the self-consistent problem for the provided region voltages.

Source code in tempura/electrostatics/solver.py

def solve(
    self,
    voltages: dict[str, float],
    initial_guess: SparseVector | None = None,
    max_ite: int = 15,
) -> None:
    """Solve the self-consistent problem for the provided region voltages."""
    sv_gate = self.pd.sparse_vector(val=0, name="Boundary")
    for key, value in voltages.items():
        sv_gate.extend(self.pd.sparse_vector(val=float(value), name=key))
    sv_charge = self.pd.sparse_vector(val=0, name="Boundary")
    poisson_input = {"voltage": sv_gate, "charge": sv_charge}
    if initial_guess is None:
        initial_guess = SparseVector(
            values=np.zeros(len(self.twodeg_idx)),
            indices=self.twodeg_idx,
            dtype=float,
        )

    self.solver.solve(
        poisson_input=poisson_input, initial_guess=initial_guess, max_ite=max_ite
    )

Convert daz (density at zero energy) to the ILDOS threshold delta.

Source code in tempura/electrostatics/solver.py

def daz_to_delta(daz: float, me_eff: float = 0.022) -> float:
    """Convert ``daz`` (density at zero energy) to the ILDOS threshold ``delta``."""
    daz = float(daz)
    return float(-daz / _helmholtz_density(me_eff))

Return ILDOS lookup tables and site labels for a 2DEG region.

The returned list contains one piecewise-linear ILDOS table per unique mesh volume in the selected 2DEG points. sites_ildos maps each 2DEG site to the corresponding ILDOS table index.

Source code in tempura/electrostatics/solver.py

def twodeg_ILDOS(
    pd: Any,
    thickness: float,
    twodeg_vol: np.ndarray,
    delta: float = 0.05,
    me_eff: float = 0.022,
    twodeg_idx: np.ndarray | None = None,
) -> tuple[list[np.ndarray], SparseVector]:
    """Return ILDOS lookup tables and site labels for a 2DEG region.

    The returned list contains one piecewise-linear ILDOS table per unique mesh
    volume in the selected 2DEG points. ``sites_ildos`` maps each 2DEG site to
    the corresponding ILDOS table index.
    """
    thickness = _require_positive_scalar(thickness, name="thickness")
    dens_helm = _helmholtz_density(me_eff)

    coord = np.empty((3, 2), dtype=float)
    coord[:, 0] = np.array([delta - 1, delta, 1])
    coord[:, 1] = np.array([0, 0, (dens_helm * (1 - delta))])

    if twodeg_idx is None:
        twodeg_idx = pd.points(name="2deg")
    twodeg_idx = np.asarray(twodeg_idx, dtype=int)
    twodeg_vol = np.asarray(twodeg_vol, dtype=float)
    if twodeg_vol.ndim != 1:
        raise ValueError(f"twodeg_vol must be 1-D; got shape {twodeg_vol.shape}.")
    if twodeg_idx.shape[0] != twodeg_vol.shape[0]:
        raise ValueError(
            "twodeg_idx and twodeg_vol must have the same length; "
            f"got {twodeg_idx.shape[0]} and {twodeg_vol.shape[0]}."
        )

    argsort = np.argsort(twodeg_vol)
    vols, counts = np.unique(np.round(twodeg_vol[argsort], 10), return_counts=True)

    sites_ildos = SparseVector(
        values=np.repeat(np.arange(len(vols), dtype=int), counts),
        indices=twodeg_idx[argsort],
    )
    ildos_arr_list: list[np.ndarray] = []
    for vol in vols:
        coord_for_volume = coord.copy()
        coord_for_volume[:, 1] *= vol / thickness
        ildos_arr_list.append(coord_for_volume)

    return ildos_arr_list, sites_ildos

Layout

Layout readers and AOI extraction helpers.

crop_polygons_to_aoi(polygons_by_layer, aoi_bbox, precision=1e-06)

Clip all layers to AOI and return local-coordinate polygons.

Source code in tempura/layout/extraction.py

def crop_polygons_to_aoi(
    polygons_by_layer: PolygonMap,
    aoi_bbox: BBoxXY,
    precision: float = 1e-6,
) -> PolygonMap:
    """Clip all layers to AOI and return local-coordinate polygons."""
    xmin, xmax, ymin, ymax = _validate_bbox(aoi_bbox, name="AOI bbox")
    clip_rect = gdstk.rectangle((xmin, ymin), (xmax, ymax))

    cropped: PolygonMap = {}
    for layer_name, polygons in polygons_by_layer.items():
        clipped = gdstk.boolean(polygons, clip_rect, "and", precision=precision)
        if not clipped:
            continue
        local_polygons: list[gdstk.Polygon] = []
        for polygon in clipped:
            points = np.asarray(polygon.points, dtype=float).copy()
            points[:, 0] -= xmin
            points[:, 1] -= ymin
            local_polygons.append(
                gdstk.Polygon(points, layer=polygon.layer, datatype=polygon.datatype)
            )
        cropped[layer_name] = local_polygons
    return cropped

load_layout_polygons(path)

Load polygons by layer from GDS/GDSII/OAS/DXF.

Source code in tempura/layout/readers.py

def load_layout_polygons(path: Path) -> PolygonMap:
    """Load polygons by layer from GDS/GDSII/OAS/DXF."""
    return load_layout_data(path)["polygons_by_layer"]

make_aoi_bbox_from_ranges(x_range, y_range)

Build AOI bbox from two [min, max]-like ranges.

Source code in tempura/layout/extraction.py

def make_aoi_bbox_from_ranges(
    x_range: Sequence[float], y_range: Sequence[float]
) -> BBoxXY:
    """Build AOI bbox from two [min, max]-like ranges."""
    if len(x_range) != 2 or len(y_range) != 2:
        raise ValueError("x_range and y_range must each have exactly 2 values.")
    xmin, xmax = sorted((float(x_range[0]), float(x_range[1])))
    ymin, ymax = sorted((float(y_range[0]), float(y_range[1])))
    return _validate_bbox((xmin, xmax, ymin, ymax), name="AOI bbox")

plot_gate_stencil_layers(gate_stencils, *, roi_size, layer_order=None, title='Rasterized gate masks by layer', figsize=None, coordinate_scale=1.0, axis_unit=None, layer_thicknesses=None)

Plot one panel per source layer with all of that layer's gates overlaid.

Source code in tempura/layout/extraction.py

def plot_gate_stencil_layers(
    gate_stencils: MasksByLayer,
    *,
    roi_size: Sequence[float],
    layer_order: Sequence[str] | None = None,
    title: str = "Rasterized gate masks by layer",
    figsize: tuple[float, float] | None = None,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
    layer_thicknesses: Mapping[str, float] | None = None,
):
    """Plot one panel per source layer with all of that layer's gates overlaid."""
    import matplotlib.pyplot as plt

    roi_length = float(roi_size[0]) * float(coordinate_scale)
    roi_width = float(roi_size[1]) * float(coordinate_scale)
    ordered_layers = sorted(gate_stencils) if layer_order is None else list(layer_order)
    non_empty_by_layer: list[tuple[str, list[np.ndarray]]] = []
    for layer_name in ordered_layers:
        masks = [
            np.asarray(mask, dtype=bool)
            for mask in gate_stencils.get(layer_name, [])
            if np.asarray(mask).any()
        ]
        if masks:
            non_empty_by_layer.append((layer_name, masks))
    if not non_empty_by_layer:
        raise ValueError("No non-empty stencils were available to plot.")

    ncols = len(non_empty_by_layer)
    if figsize is None:
        figsize = (3.6 * ncols, 3.6)
    fig, axes = plt.subplots(
        1,
        ncols,
        figsize=figsize,
        sharex=True,
        sharey=True,
    )
    axes = np.atleast_1d(axes)
    color_cycle = plt.rcParams["axes.prop_cycle"].by_key().get("color", [])

    for idx, (ax, (layer_name, masks)) in enumerate(zip(axes, non_empty_by_layer, strict=False)):
        union_mask = np.zeros_like(masks[0], dtype=bool)
        for mask in masks:
            union_mask |= mask
        ax.imshow(
            union_mask,
            origin="lower",
            cmap="Greys",
            interpolation="nearest",
            vmin=0,
            vmax=1,
            extent=(0.0, roi_length, 0.0, roi_width),
        )

        if union_mask.shape[0] > 1 and union_mask.shape[1] > 1:
            x_step = roi_length / union_mask.shape[1]
            y_step = roi_width / union_mask.shape[0]
            x_values = (np.arange(union_mask.shape[1], dtype=float) + 0.5) * x_step
            y_values = (np.arange(union_mask.shape[0], dtype=float) + 0.5) * y_step
            outline_color = (
                color_cycle[idx % len(color_cycle)] if color_cycle else "#1f77b4"
            )
            for mask in masks:
                ax.contour(
                    x_values,
                    y_values,
                    np.asarray(mask, dtype=float),
                    levels=[0.5],
                    colors=[outline_color],
                    linewidths=0.8,
                    alpha=0.85,
                )

        gate_count = len(masks)
        gate_label = "gate" if gate_count == 1 else "gates"
        title_parts = [layer_name]
        if layer_thicknesses is not None and layer_name in layer_thicknesses:
            thickness = float(layer_thicknesses[layer_name]) * float(coordinate_scale)
            unit_suffix = f" {axis_unit}" if axis_unit else ""
            title_parts.append(f"t = {thickness:.0f}{unit_suffix}")
        title_parts.append(f"{gate_count} {gate_label}")
        ax.set_title("\n".join(title_parts))
        x_label = "x" if axis_unit is None else f"x ({axis_unit})"
        y_label = "y" if axis_unit is None else f"y ({axis_unit})"
        ax.set_xlabel(x_label)
        if idx == 0:
            ax.set_ylabel(y_label)
        ax.set_aspect("equal", adjustable="box")
        ax.grid(False)

    fig.suptitle(title)
    fig.tight_layout(rect=(0.0, 0.0, 1.0, 0.95))
    return fig

plot_layout_interactive(polygons_by_layer, layers=None, exclude_layers=None, title='Layout Layers (interactive)', *, filled=True)

Return interactive Plotly figure for layer inspection.

Source code in tempura/layout/extraction.py

def plot_layout_interactive(
    polygons_by_layer: PolygonMap,
    layers: Sequence[str] | None = None,
    exclude_layers: Sequence[str] | None = None,
    title: str = "Layout Layers (interactive)",
    *,
    filled: bool = True,
) -> "go.Figure":
    """Return interactive Plotly figure for layer inspection."""
    import matplotlib.colors as mcolors
    import matplotlib.pyplot as plt
    import plotly.graph_objects as go

    layer_order = sorted(polygons_by_layer) if layers is None else list(layers)
    if exclude_layers is not None:
        excluded = {str(layer_name) for layer_name in exclude_layers}
        layer_order = [layer_name for layer_name in layer_order if layer_name not in excluded]
    fig = go.Figure()
    color_cycle = plt.rcParams["axes.prop_cycle"].by_key().get("color", [])

    for idx, layer_name in enumerate(layer_order):
        if layer_name not in polygons_by_layer:
            continue
        base_color = color_cycle[idx % len(color_cycle)] if color_cycle else "#1f77b4"
        rgba = mcolors.to_rgba(base_color, alpha=0.4)
        fill_rgba = (
            f"rgba({int(rgba[0] * 255)}, {int(rgba[1] * 255)}, "
            f"{int(rgba[2] * 255)}, {rgba[3]})"
        )
        line_rgba = (
            f"rgba({int(rgba[0] * 255)}, {int(rgba[1] * 255)}, "
            f"{int(rgba[2] * 255)}, 0.9)"
        )
        show_legend = True
        for polygon in polygons_by_layer[layer_name]:
            points = np.asarray(polygon.points, dtype=float)
            if len(points) == 0:
                continue
            if not np.allclose(points[0], points[-1]):
                points = np.vstack([points, points[0]])
            fig.add_trace(
                go.Scatter(
                    x=points[:, 0].tolist(),
                    y=points[:, 1].tolist(),
                    mode="lines",
                    fill="toself" if filled else None,
                    name=layer_name,
                    legendgroup=layer_name,
                    showlegend=show_legend,
                    line={"color": line_rgba, "width": 1},
                    fillcolor=fill_rgba if filled else "rgba(0, 0, 0, 0)",
                )
            )
            show_legend = False

    fig.update_layout(
        title=title,
        legend_title_text="Layers",
        legend={"groupclick": "togglegroup"},
        template="plotly_white",
        height=760,
        autosize=True,
    )
    fig.update_yaxes(scaleanchor="x", scaleratio=1)
    return fig

plot_layout_layers(polygons_by_layer, *, layers=None, exclude_layers=None, title=None, filled=True, coordinate_scale=1.0, axis_unit=None)

Plot layout polygons with Matplotlib using one shared set of axes.

Source code in tempura/layout/extraction.py

def plot_layout_layers(
    polygons_by_layer: PolygonMap,
    *,
    layers: Sequence[str] | None = None,
    exclude_layers: Sequence[str] | None = None,
    title: str | None = None,
    filled: bool = True,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
):
    """Plot layout polygons with Matplotlib using one shared set of axes."""
    import matplotlib.pyplot as plt

    selected_layers = sorted(polygons_by_layer)
    if layers is not None:
        allowed = set(layers)
        selected_layers = [layer for layer in selected_layers if layer in allowed]
    if exclude_layers is not None:
        excluded = set(exclude_layers)
        selected_layers = [layer for layer in selected_layers if layer not in excluded]
    if not selected_layers:
        raise ValueError("No layout layers were available to plot.")

    figure, ax = plt.subplots(figsize=(7.2, 5.8))
    color_cycle = plt.rcParams["axes.prop_cycle"].by_key().get("color", [])

    for layer_idx, layer_name in enumerate(selected_layers):
        layer_color = color_cycle[layer_idx % len(color_cycle)] if color_cycle else "#1f77b4"
        label_added = False
        for polygon in polygons_by_layer.get(layer_name, []):
            points = np.asarray(polygon.points, dtype=float)
            if len(points) == 0:
                continue
            points = points.copy()
            points[:, :2] *= float(coordinate_scale)
            if not np.allclose(points[0], points[-1]):
                points = np.vstack([points, points[0]])
            if filled:
                ax.fill(
                    points[:, 0],
                    points[:, 1],
                    linewidth=0.8,
                    alpha=0.35,
                    color=layer_color,
                    label=layer_name if not label_added else None,
                )
            else:
                ax.plot(
                    points[:, 0],
                    points[:, 1],
                    linewidth=0.9,
                    color=layer_color,
                    label=layer_name if not label_added else None,
                )
            label_added = True

    x_label = "x" if axis_unit is None else f"x ({axis_unit})"
    y_label = "y" if axis_unit is None else f"y ({axis_unit})"
    ax.set_xlabel(x_label)
    ax.set_ylabel(y_label)
    ax.set_aspect("equal", adjustable="box")
    ax.grid(False)
    if title is not None:
        ax.set_title(title)
    if selected_layers:
        ax.legend(loc="upper right", frameon=False, fontsize="small")
    figure.tight_layout()
    return figure

polygons_to_vertices(polygons_by_layer)

Convert polygons grouped by layer into nested Python lists.

Source code in tempura/layout/extraction.py

def polygons_to_vertices(polygons_by_layer: PolygonMap) -> VerticesByLayer:
    """Convert polygons grouped by layer into nested Python lists."""
    return {
        layer_name: [polygon.points.tolist() for polygon in polygons]
        for layer_name, polygons in polygons_by_layer.items()
    }

rasterize_gate_vertices(gate_vertices_by_layer, *, dx, dy, roi_bbox)

Rasterize gate polygons into boolean masks on a regular grid. Converts vector polygon data (gate vertices) organized by layer into rasterized boolean masks suitable for mask-based analysis or image processing operations. All polygons are rasterized on the same ROI grid. Args: gate_vertices_by_layer: Dictionary mapping layer names to lists of polygon vertices. Each polygon is defined by an ordered sequence of (x, y) coordinates forming a closed shape. dx: Grid spacing in the x-direction (horizontal). Must be positive. Determines the horizontal resolution of the rasterized output. dy: Grid spacing in the y-direction (vertical). Must be positive. Determines the vertical resolution of the rasterized output. roi_bbox: Rasterization region as (xmin, xmax, ymin, ymax). Every polygon mask is generated on this ROI-sized grid so all masks have identical shape and alignment. Returns: Dictionary mapping layer names to lists of boolean numpy arrays (masks). Each mask has shape (ny, nx) where a True value indicates the grid point is inside the corresponding polygon. Raises: ValueError: If dx/dy are invalid or roi_bbox bounds are invalid. Notes: - Polygons with fewer than 3 vertices are skipped (degenerate cases). - Point-in-polygon testing uses matplotlib's Path.contains_points() method. - Grid points are sampled at cell centers (offset by 0.5 * dx, 0.5 * dy). - Output masks are indexed as [y, x] (row-major ordering). - Layers with no valid polygons are omitted from the output dictionary.

Source code in tempura/layout/extraction.py

def rasterize_gate_vertices(
    gate_vertices_by_layer: VerticesByLayer,
    *,
    dx: float,
    dy: float,
    roi_bbox: BBoxXY,
) -> MasksByLayer:
    """
    Rasterize gate polygons into boolean masks on a regular grid.
    Converts vector polygon data (gate vertices) organized by layer into rasterized
    boolean masks suitable for mask-based analysis or image processing operations.
    All polygons are rasterized on the same ROI grid.
    Args:
        gate_vertices_by_layer: Dictionary mapping layer names to lists of polygon
            vertices. Each polygon is defined by an ordered sequence of (x, y) coordinates
            forming a closed shape.
        dx: Grid spacing in the x-direction (horizontal). Must be positive. Determines
            the horizontal resolution of the rasterized output.
        dy: Grid spacing in the y-direction (vertical). Must be positive. Determines
            the vertical resolution of the rasterized output.
        roi_bbox: Rasterization region as ``(xmin, xmax, ymin, ymax)``.
            Every polygon mask is generated on this ROI-sized grid so all masks have
            identical shape and alignment.
    Returns:
        Dictionary mapping layer names to lists of boolean numpy arrays (masks).
        Each mask has shape (ny, nx) where a True value indicates the grid point
        is inside the corresponding polygon.
    Raises:
        ValueError: If dx/dy are invalid or roi_bbox bounds are invalid.
    Notes:
        - Polygons with fewer than 3 vertices are skipped (degenerate cases).
        - Point-in-polygon testing uses matplotlib's Path.contains_points() method.
        - Grid points are sampled at cell centers (offset by 0.5 * dx, 0.5 * dy).
        - Output masks are indexed as [y, x] (row-major ordering).
        - Layers with no valid polygons are omitted from the output dictionary.
    """

    from matplotlib.path import Path as MplPath

    # Validate grid spacing early to avoid divide-by-zero and invalid grids.
    if dx <= 0 or dy <= 0:
        raise ValueError("dx and dy must be positive.")

    roi_xmin, roi_xmax, roi_ymin, roi_ymax = _validate_bbox(
        roi_bbox, name="roi_bbox"
    )
    roi_nx = _grid_count_from_span(roi_xmax - roi_xmin, dx)
    roi_ny = _grid_count_from_span(roi_ymax - roi_ymin, dy)
    # Sample at cell centers so the rasterized masks align with the device-grid
    # convention used throughout the electrostatics pipeline.
    roi_xs = roi_xmin + (np.arange(roi_nx) + 0.5) * dx
    roi_ys = roi_ymin + (np.arange(roi_ny) + 0.5) * dy
    roi_grid_x, roi_grid_y = np.meshgrid(roi_xs, roi_ys)
    roi_sample_points = np.column_stack([roi_grid_x.ravel(), roi_grid_y.ravel()])

    # Accumulate per-layer masks; omit layers with no valid polygons.
    masks_by_layer: MasksByLayer = {}
    for layer_name, polygons in gate_vertices_by_layer.items():
        # Collect masks for each polygon in the current layer.
        layer_masks: list[np.ndarray] = []
        for vertices in polygons:
            # Normalize vertex data and skip degenerate polygons.
            points = np.asarray(vertices, dtype=float)
            if len(points) < 3:
                continue
            # Rasterize on the shared ROI grid to preserve global alignment.
            mask = MplPath(points).contains_points(roi_sample_points).reshape(
                roi_ny, roi_nx
            )
            # Append the rasterized mask for this polygon.
            layer_masks.append(mask)
        if layer_masks:
            # Keep only layers that have at least one valid mask.
            masks_by_layer[layer_name] = layer_masks
    return masks_by_layer

read_dxf_polylines(path)

Read closed DXF POLYLINE entities grouped by DXF layer name.

Source code in tempura/layout/readers.py

def read_dxf_polylines(path: Path) -> PolygonMap:
    """Read closed DXF ``POLYLINE`` entities grouped by DXF layer name."""
    return read_dxf_layout(path)["polygons_by_layer"]

read_gds_polygons(path)

Read polygons from a GDS/GDSII layout.

Source code in tempura/layout/readers.py

def read_gds_polygons(path: Path) -> PolygonMap:
    """Read polygons from a GDS/GDSII layout."""
    return read_gds_layout(path)["polygons_by_layer"]

read_oas_polygons(path)

Read polygons from an OASIS layout.

Source code in tempura/layout/readers.py

def read_oas_polygons(path: Path) -> PolygonMap:
    """Read polygons from an OASIS layout."""
    return read_oas_layout(path)["polygons_by_layer"]

Load polygons and available unit metadata from GDS/GDSII/OAS/DXF.

Source code in tempura/layout/readers.py

def load_layout_data(path: Path) -> LayoutData:
    """Load polygons and available unit metadata from GDS/GDSII/OAS/DXF."""
    if not path.exists():
        raise FileNotFoundError(f"Layout file not found: {path}")

    suffix = path.suffix.lower()
    try:
        reader = LAYOUT_READERS[suffix]
    except KeyError as exc:
        raise ValueError(
            f"Unsupported layout suffix {suffix!r}; use {sorted(SUPPORTED_LAYOUT_SUFFIXES)}"
        ) from exc
    return reader(path)

Load a layout file and return an interactive figure for its polygons.

Source code in tempura/layout/extraction.py

def plot_layout_file_interactive(
    layout_path: str | Path,
    *,
    layers: Sequence[str] | None = None,
    exclude_layers: Sequence[str] | None = None,
    title: str | None = None,
    filled: bool = True,
) -> "go.Figure":
    """Load a layout file and return an interactive figure for its polygons."""
    from tempura.layout.readers import load_layout_data

    path = Path(layout_path)
    layout_data = load_layout_data(path)
    if title is None:
        title = f"{path.name} layout"
    return plot_layout_interactive(
        layout_data["polygons_by_layer"],
        layers=layers,
        exclude_layers=exclude_layers,
        title=title,
        filled=filled,
    )

Shared helpers for layout-backed demo simulations.

build_device(prepared, stack)

Build a finalized device from prepared masks and an ordered stack spec.

The stack describes the physical layers and may optionally include one per-layer resolution entry used directly by build_problem(...). 2deg entries may also provide boundary_condition to select the linear response model for that region. Gate layers may also set inverted=True when their one source stencil represents an etched opening rather than deposited metal.

Source code in tempura/layout/layout_pipeline.py

def build_device(
    prepared: dict[str, object],
    stack: list[dict[str, object]],
) -> Device:
    """Build a finalized device from prepared masks and an ordered stack spec.

    The stack describes the physical layers and may optionally include one
    per-layer ``resolution`` entry used directly by ``build_problem(...)``.
    ``2deg`` entries may also provide ``boundary_condition`` to select the
    linear response model for that region. Gate layers may also set
    ``inverted=True`` when their one source stencil represents an etched
    opening rather than deposited metal.
    """
    grid_step = _require_positive(float(prepared["grid_step"]), name="grid_step")
    roi_length, roi_width = prepared["roi_size"]
    roi_length = _require_positive(float(roi_length), name="roi_length")
    roi_width = _require_positive(float(roi_width), name="roi_width")
    expected_shape = tuple(int(v) for v in prepared["grid_shape"])
    gate_stencils = prepared["gate_stencils"]

    device = Device(length=roi_length, width=roi_width, grid=grid_step)
    for layer in stack:
        kind = str(layer["kind"])
        name = str(layer["name"])
        if "invert" in layer:
            raise ValueError(f"{name} uses 'invert'; use 'inverted' instead.")
        inverted = bool(layer.get("inverted", False))
        source_layer_value = layer.get("source_layer")
        if inverted and source_layer_value is None:
            raise ValueError(
                f"{name}.inverted requires source_layer because there is "
                "no stencil to invert."
            )
        if inverted and kind != "gate":
            raise ValueError(f"{name}.inverted is only supported for gate layers.")
        thickness = _require_integer_size(
            float(layer["thickness"]),
            name=f"{name}.thickness",
        )
        resolution = _layer_resolution_from_spec(layer, grid_step=grid_step)

        if kind == "dielectric":
            device.add_layer(
                Dielectric(
                    name=name,
                    permittivity=float(layer["permittivity"]),
                    thickness=thickness,
                    resolution=resolution,
                )
            )
            continue

        if kind == "2deg":
            device.add_layer(
                TwoDEG(
                    name=name,
                    permittivity=float(layer["permittivity"]),
                    thickness=thickness,
                    resolution=resolution,
                    boundary_condition=str(layer.get("boundary_condition", "helmholtz")),
                )
            )
            continue

        if kind != "gate":
            raise ValueError(f"Unsupported stack layer kind {kind!r}.")

        if source_layer_value is None:
            raise ValueError(f"{name} requires source_layer.")
        source_layer = str(source_layer_value)
        masks = gate_stencils.get(source_layer)
        if not masks:
            raise ValueError(f"Missing gate stencils for source_layer {source_layer!r}.")
        source_masks: list[np.ndarray] = []
        for gate_idx, stencil in enumerate(masks):
            mask = np.asarray(stencil, dtype=bool)
            if mask.shape != expected_shape:
                raise ValueError(
                    f"Stencil shape mismatch for {source_layer}[{gate_idx}]: "
                    f"expected {expected_shape}, got {mask.shape}."
                )
            if np.any(mask):
                source_masks.append(mask)

        if inverted:
            if len(source_masks) != 1:
                raise ValueError(
                    f"{name}.inverted requires exactly one non-empty stencil in "
                    f"source_layer {source_layer!r}; got {len(source_masks)}."
                )
            component_masks = _split_connected_stencils(~source_masks[0])
        else:
            component_masks = source_masks

        gate_batch: list[Gate] = []
        for mask in component_masks:
            # Layout-backed gate layers expand into one Device gate per connected
            # rasterized region. In inverted mode those regions come from the
            # connected components of one inverted opening mask.
            # Cropping and rasterization can leave behind empty connected
            # components at the AOI boundary. They should not become Device gates.
            if not np.any(mask):
                continue
            gate_batch.append(
                Gate(
                    name=f"{name}_{len(gate_batch)}",
                    thickness=thickness,
                    stencil=mask,
                    resolution=resolution,
                )
            )

        if not gate_batch:
            raise ValueError(
                f"Source layer {source_layer!r} did not produce any non-empty gate masks."
            )
        device.add_layer(gate_batch)

    device.finalize()
    return device

format_device_dimensions(device, prepared=None)

Return the physical size and XY grid summary for a device.

Source code in tempura/layout/layout_pipeline.py

def format_device_dimensions(
    device: Device,
    prepared: dict[str, object] | None = None,
) -> str:
    """Return the physical size and XY grid summary for a device."""
    ny, nx = tuple(int(v) for v in device.grid_shape)
    summary = {
        "Lx": device.length,
        "Ly": device.width,
        "grid": device.grid,
        "grid_shape": f"(ny={ny}, nx={nx})",
        "grid_points": ny * nx,
        "layers": ", ".join(device.layers),
    }
    if prepared is not None:
        roi_size_m = prepared.get("roi_size_m")
        grid_constant_m = prepared.get("grid_constant_m")
        unit_scale_m = prepared.get("layout_unit_scale_m")
        if roi_size_m is not None:
            summary["physical_size"] = (
                f"({_format_meters(roi_size_m[0])}, {_format_meters(roi_size_m[1])})"
            )
        if grid_constant_m is not None:
            summary["grid_constant"] = _format_meters(grid_constant_m)
        if unit_scale_m is not None:
            summary["layout_unit_scale"] = _format_meters(unit_scale_m)
    return format_mapping_block("Device summary", summary)

format_roi_summary(prepared)

Return the ROI size in physical units and grid-normalized units.

Source code in tempura/layout/layout_pipeline.py

def format_roi_summary(prepared: dict[str, object]) -> str:
    """Return the ROI size in physical units and grid-normalized units."""
    roi_length_m, roi_width_m = prepared["roi_size_m"]
    roi_length, roi_width = prepared["roi_size"]
    ny, nx = tuple(int(v) for v in prepared["grid_shape"])
    grid_constant_m = float(prepared["grid_constant_m"])
    layout_width, layout_height = prepared["roi_size_layout"]

    summary = {
        "layout_size": f"({layout_width}, {layout_height})",
        "physical_size": f"({_format_meters(roi_length_m)}, {_format_meters(roi_width_m)})",
        "grid_constant": _format_meters(grid_constant_m),
        "rescaled_size": f"({roi_length}, {roi_width})",
        "grid_shape": f"(ny={ny}, nx={nx})",
        "grid_points": ny * nx,
    }
    unit_scale_m = prepared.get("layout_unit_scale_m")
    if unit_scale_m is not None:
        summary["layout_unit_scale"] = _format_meters(unit_scale_m)
    return format_mapping_block("ROI summary", summary)

prepare_layout(layout_path, aoi_bbox, size_mode, grid_constant_m=DEFAULT_GRID_CONSTANT_M, physical_x_length=None, precision=1e-06, tolerance=1e-06)

Load, crop, scale, and rasterize a layout-backed ROI.

The ROI is first resolved in physical units (meters), then rescaled so one internal XY unit corresponds to one grid constant.

Source code in tempura/layout/layout_pipeline.py

def prepare_layout(
    layout_path: str | Path,
    aoi_bbox: tuple[float, float, float, float],
    size_mode: str,
    grid_constant_m: float = DEFAULT_GRID_CONSTANT_M,
    physical_x_length: float | None = None,
    precision: float = 1e-6,
    tolerance: float = 1e-6,
) -> dict[str, object]:
    """Load, crop, scale, and rasterize a layout-backed ROI.

    The ROI is first resolved in physical units (meters), then rescaled so one
    internal XY unit corresponds to one grid constant.
    """
    grid_constant_m = _require_positive(grid_constant_m, name="grid_constant_m")

    layout_data = load_layout_data(Path(layout_path))
    polygons_by_layer = layout_data["polygons_by_layer"]
    cropped_polygons_by_layer = crop_polygons_to_aoi(
        polygons_by_layer,
        aoi_bbox,
        precision=precision,
    )

    layout_width = float(aoi_bbox[1] - aoi_bbox[0])
    layout_height = float(aoi_bbox[3] - aoi_bbox[2])
    if layout_width <= 0 or layout_height <= 0:
        raise ValueError("aoi_bbox must define a positive-width, positive-height ROI.")

    if size_mode == "layout_size":
        if not bool(layout_data.get("units_known", False)):
            raise ValueError(
                "size_mode='layout_size' requires layout unit metadata, but none "
                f"was found for {layout_path}."
            )
        unit_scale_m = float(layout_data["unit_scale_m"])
        roi_length_m = layout_width * unit_scale_m
        roi_width_m = layout_height * unit_scale_m
    elif size_mode == "explicit_x":
        if physical_x_length is None:
            raise ValueError("physical_x_length is required for size_mode='explicit_x'.")
        roi_length_m = _require_positive(physical_x_length, name="physical_x_length")
        roi_width_m = roi_length_m * (layout_height / layout_width)
    else:
        raise ValueError(
            "size_mode must be either 'layout_size' or 'explicit_x'; "
            f"got {size_mode!r}."
        )

    roi_length_internal, nx = _require_grid_multiple(
        roi_length_m,
        grid_constant_m,
        name="roi_length",
        tolerance=tolerance,
    )
    roi_width_internal, ny = _require_grid_multiple(
        roi_width_m,
        grid_constant_m,
        name="roi_width",
        tolerance=tolerance,
    )

    scaled_cropped_polygons = _scale_polygons(
        cropped_polygons_by_layer,
        sx=roi_length_internal / layout_width,
        sy=roi_width_internal / layout_height,
    )
    gate_vertices_by_layer = polygons_to_vertices(scaled_cropped_polygons)
    gate_stencils = rasterize_gate_vertices(
        gate_vertices_by_layer,
        dx=1.0,
        dy=1.0,
        roi_bbox=(0.0, roi_length_internal, 0.0, roi_width_internal),
    )
    unit_scale_m = layout_data.get("unit_scale_m")
    roi_size_m = (roi_length_m, roi_width_m)

    return {
        "polygons_by_layer": polygons_by_layer,
        "cropped_polygons_by_layer": scaled_cropped_polygons,
        "gate_vertices_by_layer": gate_vertices_by_layer,
        "gate_stencils": gate_stencils,
        "roi_size": (roi_length_internal, roi_width_internal),
        "grid_shape": (ny, nx),
        "grid_step": 1.0,
        "grid_constant_m": grid_constant_m,
        "aoi_bbox": aoi_bbox,
        "layout_unit_scale_m": unit_scale_m,
        "layout_unit_source": layout_data.get("unit_source"),
        "layout_units_known": bool(layout_data.get("units_known", False)),
        "roi_size_m": roi_size_m,
        "roi_size_layout": (layout_width, layout_height),
    }

print_device_dimensions(device, prepared=None)

Print a compact device summary.

Source code in tempura/layout/layout_pipeline.py

def print_device_dimensions(
    device: Device,
    prepared: dict[str, object] | None = None,
) -> None:
    """Print a compact device summary."""
    print(format_device_dimensions(device, prepared))

print_roi_summary(prepared)

Print a compact ROI summary.

Source code in tempura/layout/layout_pipeline.py

def print_roi_summary(prepared: dict[str, object]) -> None:
    """Print a compact ROI summary."""
    print(format_roi_summary(prepared))

solve_demo(device, vacuum_scale=2, rhs_block_size=1)

Solve the electrostatic basis problem and extract a 2DEG slice.

Source code in tempura/layout/layout_pipeline.py

def solve_demo(
    device: Device,
    vacuum_scale: float = 2,
    rhs_block_size: int = 1,
) -> dict[str, object]:
    """Solve the electrostatic basis problem and extract a 2DEG slice."""
    pb, region_shapes, gate_names = build_problem(
        device,
        vacuum_scale=vacuum_scale,
    )
    basis_potentials, problem = solve_gate_potentials(
        pb,
        region_shapes,
        gate_names,
        rhs_block_size=rhs_block_size,
    )

    plane = extract_2deg_plane(problem, region_shapes)
    if (plane.ny, plane.nx) != tuple(int(v) for v in device.grid_shape):
        raise ValueError(
            "2DEG grid shape mismatch: "
            f"expected {device.grid_shape}, got {(plane.ny, plane.nx)}."
        )

    basis_2deg = extract_2deg_basis(basis_potentials, gate_names, plane)

    return {
        "gate_names": gate_names,
        "basis_potentials": basis_potentials,
        "basis_2deg": basis_2deg,
        "coords_2deg": plane.coords,
        "problem": problem,
        "region_shapes": region_shapes,
    }

Plotting

Plotting helpers for Tempura.

PlaneSpec dataclass

Describe one 2D slice through a finalized 3D problem.

Attributes:

Name	Type	Description
title	str	Panel title.
plot_region	Box	2D plotting region in the displayed axes.
plane_axis	int	Fixed axis in 3D coordinates.
plane_value	float	Position of the fixed slicing plane.
axis_indices	tuple[int, int]	Indices of the displayed axes in the 3D coordinate system.
axis_labels	tuple[str, str]	Labels for the displayed axes.

Source code in tempura/plotting/meshing.py

@dataclass(frozen=True)
class PlaneSpec:
    """Describe one 2D slice through a finalized 3D problem.

    Attributes:
        title: Panel title.
        plot_region: 2D plotting region in the displayed axes.
        plane_axis: Fixed axis in 3D coordinates.
        plane_value: Position of the fixed slicing plane.
        axis_indices: Indices of the displayed axes in the 3D coordinate system.
        axis_labels: Labels for the displayed axes.
    """

    title: str
    plot_region: shapes.Box
    plane_axis: int
    plane_value: float
    axis_indices: tuple[int, int]
    axis_labels: tuple[str, str]

    def __post_init__(self) -> None:
        bbox = np.asarray(self.plot_region.bbox, dtype=float)
        if bbox.shape != (2, 2):
            raise ValueError(
                "PlaneSpec.plot_region must be a 2D shape with bbox shape (2, 2); "
                f"got {bbox.shape}."
            )
        if len(self.axis_indices) != 2 or len(set(self.axis_indices)) != 2:
            raise ValueError(
                "PlaneSpec.axis_indices must contain exactly two distinct axes; "
                f"got {self.axis_indices}."
            )
        if self.plane_axis in self.axis_indices:
            raise ValueError(
                "PlaneSpec.plane_axis must differ from the plotted axis indices; "
                f"got plane_axis={self.plane_axis}, axis_indices={self.axis_indices}."
            )
        if len(self.axis_labels) != 2:
            raise ValueError(
                "PlaneSpec.axis_labels must contain exactly two strings; "
                f"got {self.axis_labels}."
            )

add_gate_footprint_contours_xy(ax, x_values, y_values, footprints, *, coordinate_scale=1.0, color='black', linewidth=0.8, alpha=0.7)

Overlay direct XY gate-footprint contours on one existing axis.

Source code in tempura/plotting/meshing.py

def add_gate_footprint_contours_xy(
    ax: plt.Axes,
    x_values: np.ndarray,
    y_values: np.ndarray,
    footprints: Mapping[str, np.ndarray],
    *,
    coordinate_scale: float = 1.0,
    color: str = "black",
    linewidth: float = 0.8,
    alpha: float = 0.7,
) -> None:
    """Overlay direct XY gate-footprint contours on one existing axis."""
    if not footprints:
        return
    scale = float(coordinate_scale)
    x_vals = np.asarray(x_values, dtype=float) * scale
    y_vals = np.asarray(y_values, dtype=float) * scale
    for mask in footprints.values():
        arr = np.asarray(mask, dtype=bool)
        if not np.any(arr):
            continue
        ax.contour(
            x_vals,
            y_vals,
            arr.astype(float),
            levels=[0.5],
            colors=[color],
            linewidths=linewidth,
            alpha=alpha,
            zorder=2,
        )

add_plane_cut_guides(ax, *, x=None, y=None, coordinate_scale=1.0, color='black', linewidth=1.0, linestyle='--')

Overlay one or two cut-guide lines on an existing 2D plane plot.

Source code in tempura/plotting/meshing.py

def add_plane_cut_guides(
    ax: plt.Axes,
    *,
    x: float | None = None,
    y: float | None = None,
    coordinate_scale: float = 1.0,
    color: str = "black",
    linewidth: float = 1.0,
    linestyle: str = "--",
) -> None:
    """Overlay one or two cut-guide lines on an existing 2D plane plot."""
    scale = float(coordinate_scale)
    if x is not None:
        ax.axvline(float(x) * scale, color=color, linewidth=linewidth, linestyle=linestyle)
    if y is not None:
        ax.axhline(float(y) * scale, color=color, linewidth=linewidth, linestyle=linestyle)

add_region_contours_on_planes(problem, axes, plane_specs, region_shapes, *, region_display_names=None, coordinate_scale=1.0, linewidth=1.0, alpha=0.9)

Overlay region-outline contours on existing plane plots.

Source code in tempura/plotting/meshing.py

def add_region_contours_on_planes(
    problem: Any,
    axes: Sequence[plt.Axes],
    plane_specs: Sequence[PlaneSpec],
    region_shapes: RegionShapeMap,
    *,
    region_display_names: RegionDisplayMap | None = None,
    coordinate_scale: float = 1.0,
    linewidth: float = 1.0,
    alpha: float = 0.9,
) -> None:
    """Overlay region-outline contours on existing plane plots."""
    if not plane_specs or not region_shapes:
        return

    coordinates = np.asarray(problem.coordinates, dtype=float)
    volume_bounds = _volume_bounds_lookups(coordinates)
    display_names, _ = _display_order(list(region_shapes), region_display_names)
    grouped_names: dict[str, list[str]] = {}
    for region_name in region_shapes:
        display_name = (
            str(region_display_names[region_name])
            if region_display_names is not None and region_name in region_display_names
            else str(region_name)
        )
        grouped_names.setdefault(display_name, []).append(region_name)

    cmap = plt.cm.get_cmap("tab10")
    colors = cmap(np.linspace(0.0, 1.0, max(len(display_names), 1)))

    for ax, spec in zip(axes, plane_specs, strict=False):
        plane_sites = _plane_sites(
            coordinates,
            spec=spec,
            bounds_lookup=volume_bounds[spec.plane_axis],
        )
        if plane_sites.size == 0:
            continue
        x_edges, y_edges = _region_sample_edges(
            plane_sites,
            spec=spec,
            volume_bounds=volume_bounds,
        )
        x_centers = 0.5 * (x_edges[:-1] + x_edges[1:])
        y_centers = 0.5 * (y_edges[:-1] + y_edges[1:])
        xx, yy = np.meshgrid(x_centers, y_centers, indexing="xy")

        sample_points = np.zeros((xx.size, 3), dtype=float)
        sample_points[:, spec.plane_axis] = float(spec.plane_value)
        sample_points[:, spec.axis_indices[0]] = xx.ravel()
        sample_points[:, spec.axis_indices[1]] = yy.ravel()

        for color_index, display_name in enumerate(display_names):
            union = np.zeros(sample_points.shape[0], dtype=bool)
            for region_name in grouped_names.get(display_name, []):
                union |= np.asarray(region_shapes[region_name](sample_points), dtype=bool)
            if not np.any(union):
                continue
            ax.contour(
                xx * float(coordinate_scale),
                yy * float(coordinate_scale),
                union.reshape(len(y_centers), len(x_centers)).astype(float),
                levels=[0.5],
                colors=[colors[color_index]],
                linewidths=linewidth,
                alpha=alpha,
                zorder=2,
            )

expand_plane_grid(plane_specs_by_key, grid, *, title_overrides=None)

Expand a named subplot grid into an ordered plane-spec sequence.

Source code in tempura/plotting/meshing.py

def expand_plane_grid(
    plane_specs_by_key: Mapping[str, PlaneSpec],
    grid: Sequence[Sequence[str]],
    *,
    title_overrides: Mapping[tuple[int, int], str] | None = None,
) -> tuple[list[PlaneSpec], int]:
    """Expand a named subplot grid into an ordered plane-spec sequence."""
    ordered_specs: list[PlaneSpec] = []
    ncols = 0
    for row_idx, row in enumerate(grid):
        ncols = max(ncols, len(row))
        for col_idx, key in enumerate(row):
            try:
                spec = plane_specs_by_key[key]
            except KeyError as exc:
                raise ValueError(f"Unknown plane key {key!r} in grid.") from exc
            title = None if title_overrides is None else title_overrides.get((row_idx, col_idx))
            ordered_specs.append(spec if title is None else replace(spec, title=title))
    if not ordered_specs:
        raise ValueError("grid must contain at least one plane key.")
    return ordered_specs, ncols

make_standard_plane_specs(device, *, xy_title, xy_plane_z, xz_title, xz_plane_y, yz_title=None, yz_plane_x=None, z_padding=1.0)

Build the standard XY/XZ/YZ plotting cuts for one device.

Source code in tempura/plotting/meshing.py

def make_standard_plane_specs(
    device: Any,
    *,
    xy_title: str,
    xy_plane_z: float,
    xz_title: str,
    xz_plane_y: float,
    yz_title: str | None = None,
    yz_plane_x: float | None = None,
    z_padding: float = 1.0,
) -> dict[str, PlaneSpec]:
    """Build the standard XY/XZ/YZ plotting cuts for one device."""
    z_values = []
    for layer in device.layers.values():
        bbox = np.asarray(layer.shape.bbox, dtype=float)
        z_values.extend([float(bbox[0, 2]), float(bbox[1, 2])])
    if not z_values:
        raise ValueError("The device does not contain any finalized layers.")

    z_min = min(z_values) - float(z_padding) / 2.0
    z_size = (max(z_values) - min(z_values)) + float(z_padding)

    specs = {
        "XY": PlaneSpec(
            title=xy_title,
            plot_region=shapes.Box(
                lower_left=np.array([-device.length / 2.0, -device.width / 2.0], dtype=float),
                size=np.array([device.length, device.width], dtype=float),
            ),
            plane_axis=2,
            plane_value=float(xy_plane_z),
            axis_indices=(0, 1),
            axis_labels=("x", "y"),
        ),
        "XZ": PlaneSpec(
            title=xz_title,
            plot_region=shapes.Box(
                lower_left=np.array([-device.length / 2.0, z_min], dtype=float),
                size=np.array([device.length, z_size], dtype=float),
            ),
            plane_axis=1,
            plane_value=float(xz_plane_y),
            axis_indices=(0, 2),
            axis_labels=("x", "z"),
        ),
    }
    if yz_title is not None and yz_plane_x is not None:
        specs["YZ"] = PlaneSpec(
            title=yz_title,
            plot_region=shapes.Box(
                lower_left=np.array([-device.width / 2.0, z_min], dtype=float),
                size=np.array([device.width, z_size], dtype=float),
            ),
            plane_axis=0,
            plane_value=float(yz_plane_x),
            axis_indices=(1, 2),
            axis_labels=("y", "z"),
        )
    return specs

make_xy_emphasis_axes(*, figsize=(12.0, 9.8), layout='stack', width_ratios=(1.0, 1.0), height_ratios=(1.0, 1.0, 1.0))

Return three axes for XY, XZ, and YZ cuts in one of two compact layouts.

Source code in tempura/plotting/meshing.py

def make_xy_emphasis_axes(
    *,
    figsize: tuple[float, float] = (12.0, 9.8),
    layout: str = "stack",
    width_ratios: tuple[float, float] = (1.0, 1.0),
    height_ratios: tuple[float, float, float] = (1.0, 1.0, 1.0),
) -> tuple[plt.Figure, np.ndarray]:
    """Return three axes for XY, XZ, and YZ cuts in one of two compact layouts."""
    fig = plt.figure(figsize=figsize, constrained_layout=True)
    if layout == "stack":
        gs = fig.add_gridspec(3, 1, height_ratios=height_ratios)
        xy_ax = fig.add_subplot(gs[0, 0])
        xz_ax = fig.add_subplot(gs[1, 0])
        yz_ax = fig.add_subplot(gs[2, 0])
    elif layout == "left_span":
        gs = fig.add_gridspec(2, 2, width_ratios=width_ratios)
        xy_ax = fig.add_subplot(gs[:, 0])
        xz_ax = fig.add_subplot(gs[0, 1])
        yz_ax = fig.add_subplot(gs[1, 1])
    else:
        raise ValueError(f"Unknown layout {layout!r}; expected 'stack' or 'left_span'.")
    return fig, np.asarray([xy_ax, xz_ax, yz_ax], dtype=object)

nearest_axis_value(coordinates, axis, target)

Snap a requested cut to the nearest finalized mesh plane.

Source code in tempura/plotting/meshing.py

def nearest_axis_value(coordinates: np.ndarray, axis: int, target: float) -> float:
    """Snap a requested cut to the nearest finalized mesh plane."""
    values = np.unique(np.round(np.asarray(coordinates, dtype=float)[:, axis], 12))
    return float(values[np.argmin(np.abs(values - target))])

plot_problem_regions_with_mesh(problem, region_shapes, plane_specs, *, axes=None, region_names=None, region_display_names=None, ncols=None, figsize=(13.0, 11.0), suptitle='Pescado region cuts with overlaid mesh nodes', fill_mode='regions', background_alpha=0.88, tile_edgecolors='white', tile_linewidth=0.6, show_volume_edges=False, show_mesh_points=True, mesh_point_size=12.0, mesh_facecolors='white', mesh_edgecolors='black', mesh_linewidth=0.35, coordinate_scale=1.0, axis_unit=None)

Plot finalized region slices with overlaid mesh nodes.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem with a coordinates array.	required
region_shapes	RegionShapeMap	Mapping returned by build_problem(...). These are the solver-owned region shapes after Tempura's lower-face trimming.	required
plane_specs	Sequence[PlaneSpec]	Sequence of PlaneSpec describing the cuts to draw.	required
region_names	Sequence[str] | None	Optional region order for coloring. Defaults to every region except "Boundary" in insertion order.	None
region_display_names	RegionDisplayMap | None	Optional mapping from solver-owned region names to display labels. Regions that share a display label are drawn with the same color and collapsed into one legend entry.	None
ncols	int | None	Number of subplot columns. Defaults to 2 when more than one plane is requested and 1 otherwise.	None
figsize	tuple[float, float]	Matplotlib figure size.	(13.0, 11.0)
suptitle	str	Figure title.	'Pescado region cuts with overlaid mesh nodes'
fill_mode	str	"regions" to render the material cross-section from the finalized region shapes, or "control_volumes" to render the selected mesh control volumes directly.	'regions'
background_alpha	float	Alpha used for the ownership tiles under the markers.	0.88
tile_edgecolors	str	Edge color for the projected mesh control volumes.	'white'
tile_linewidth	float	Edge line width for the projected mesh control volumes.	0.6
show_volume_edges	bool	Whether to overlay the projected control-volume boundaries on top of the filled region cross-section.	False
show_mesh_points	bool	Whether to overlay the selected mesh sites as points.	True
mesh_point_size	float	Scatter marker size for finalized mesh nodes.	12.0
mesh_facecolors	str	Scatter face color.	'white'
mesh_edgecolors	str	Scatter edge color.	'black'
mesh_linewidth	float	Scatter edge line width.	0.35
coordinate_scale	float	Multiplicative scale applied to displayed coordinates.	1.0
axis_unit	str | None	Optional axis-unit label appended to x/y/z labels.	None

Returns:

Type	Description
tuple[Figure, ndarray]	Tuple of (figure, axes).

Notes

For each requested plane, Tempura first selects the mesh sites whose control volumes intersect that plane. The optional marker overlay shows those selected sites explicitly. The colored background can either show the projected control volumes directly or a sampled cross-section of the finalized region shapes.

Source code in tempura/plotting/meshing.py

def plot_problem_regions_with_mesh(
    problem: Any,
    region_shapes: RegionShapeMap,
    plane_specs: Sequence[PlaneSpec],
    *,
    axes: Sequence[plt.Axes] | None = None,
    region_names: Sequence[str] | None = None,
    region_display_names: RegionDisplayMap | None = None,
    ncols: int | None = None,
    figsize: tuple[float, float] = (13.0, 11.0),
    suptitle: str = "Pescado region cuts with overlaid mesh nodes",
    fill_mode: str = "regions",
    background_alpha: float = 0.88,
    tile_edgecolors: str = "white",
    tile_linewidth: float = 0.6,
    show_volume_edges: bool = False,
    show_mesh_points: bool = True,
    mesh_point_size: float = 12.0,
    mesh_facecolors: str = "white",
    mesh_edgecolors: str = "black",
    mesh_linewidth: float = 0.35,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
) -> tuple[plt.Figure, np.ndarray]:
    """Plot finalized region slices with overlaid mesh nodes.

    Args:
        problem: Finalized Pescado problem with a ``coordinates`` array.
        region_shapes: Mapping returned by ``build_problem(...)``. These are the
            solver-owned region shapes after Tempura's lower-face trimming.
        plane_specs: Sequence of ``PlaneSpec`` describing the cuts to draw.
        region_names: Optional region order for coloring. Defaults to every
            region except ``"Boundary"`` in insertion order.
        region_display_names: Optional mapping from solver-owned region names to
            display labels. Regions that share a display label are drawn with
            the same color and collapsed into one legend entry.
        ncols: Number of subplot columns. Defaults to ``2`` when more than one
            plane is requested and ``1`` otherwise.
        figsize: Matplotlib figure size.
        suptitle: Figure title.
        fill_mode: ``"regions"`` to render the material cross-section from the
            finalized region shapes, or ``"control_volumes"`` to render the
            selected mesh control volumes directly.
        background_alpha: Alpha used for the ownership tiles under the markers.
        tile_edgecolors: Edge color for the projected mesh control volumes.
        tile_linewidth: Edge line width for the projected mesh control volumes.
        show_volume_edges: Whether to overlay the projected control-volume
            boundaries on top of the filled region cross-section.
        show_mesh_points: Whether to overlay the selected mesh sites as points.
        mesh_point_size: Scatter marker size for finalized mesh nodes.
        mesh_facecolors: Scatter face color.
        mesh_edgecolors: Scatter edge color.
        mesh_linewidth: Scatter edge line width.
        coordinate_scale: Multiplicative scale applied to displayed coordinates.
        axis_unit: Optional axis-unit label appended to x/y/z labels.

    Returns:
        Tuple of ``(figure, axes)``.

    Notes:
        For each requested plane, Tempura first selects the mesh sites whose
        control volumes intersect that plane. The optional marker overlay shows
        those selected sites explicitly. The colored background can either show
        the projected control volumes directly or a sampled cross-section of
        the finalized region shapes.
    """
    if not plane_specs:
        raise ValueError("plane_specs must not be empty.")
    if fill_mode not in {"regions", "control_volumes"}:
        raise ValueError(
            "fill_mode must be either 'regions' or 'control_volumes'; "
            f"got {fill_mode!r}."
        )
    if ncols is None:
        ncols = 2 if len(plane_specs) > 1 else 1
    if ncols <= 0:
        raise ValueError(f"ncols must be positive; got {ncols}.")

    coordinates = np.asarray(problem.coordinates, dtype=float)
    if coordinates.ndim != 2 or coordinates.shape[1] != 3:
        raise ValueError(
            "problem.coordinates must have shape (npoints, 3); "
            f"got {coordinates.shape}."
        )

    if region_names is None:
        region_names = [name for name in region_shapes if name != "Boundary"]
    else:
        region_names = list(region_names)
    if not region_names:
        raise ValueError("No region names were provided for plotting.")
    display_names, display_lookup = _display_order(region_names, region_display_names)

    volume_bounds = _volume_bounds_lookups(coordinates)
    if axes is None:
        nrows = int(np.ceil(len(plane_specs) / ncols))
        fig, axes_obj = plt.subplots(nrows, ncols, figsize=figsize, constrained_layout=True)
        axes_array = np.atleast_1d(axes_obj).ravel()
        created_axes = True
    else:
        axes_array = np.atleast_1d(np.asarray(list(axes), dtype=object)).ravel()
        if len(axes_array) < len(plane_specs):
            raise ValueError(
                f"axes must contain at least {len(plane_specs)} axes; got {len(axes_array)}."
            )
        fig = axes_array[0].figure
        created_axes = False

    for ax, spec in zip(axes_array, plane_specs, strict=False):
        plot_bbox = np.asarray(spec.plot_region.bbox, dtype=float) * float(coordinate_scale)
        plane_sites = _plane_sites(
            coordinates,
            spec=spec,
            bounds_lookup=volume_bounds[spec.plane_axis],
        )
        plane_points = (
            np.asarray(plane_sites[:, spec.axis_indices], dtype=float) * float(coordinate_scale)
        )
        labels = _label_plane_sites(
            plane_sites,
            region_shapes,
            region_names,
            region_display_names=region_display_names,
            display_lookup=display_lookup,
        )
        if fill_mode == "regions":
            _region_fill(
                ax,
                spec=spec,
                region_shapes=region_shapes,
                region_names=region_names,
                region_display_names=region_display_names,
                display_lookup=display_lookup,
                display_names=display_names,
                plane_sites=plane_sites,
                volume_bounds=volume_bounds,
                background_alpha=background_alpha,
                coordinate_scale=coordinate_scale,
            )
            if show_volume_edges:
                volume_edges = _volume_edges(
                    plane_sites,
                    spec=spec,
                    coordinate_scale=coordinate_scale,
                    volume_bounds=volume_bounds,
                    tile_edgecolors=tile_edgecolors,
                    tile_linewidth=tile_linewidth,
                )
                if volume_edges is not None:
                    ax.add_collection(volume_edges)
        else:
            ownership_tiles = _ownership_tiles(
                plane_sites,
                labels,
                spec=spec,
                coordinate_scale=coordinate_scale,
                background_alpha=background_alpha,
                volume_bounds=volume_bounds,
                tile_edgecolors=tile_edgecolors,
                tile_linewidth=tile_linewidth,
            )
            if ownership_tiles is not None:
                ax.add_collection(ownership_tiles)
        if show_mesh_points:
            ax.scatter(
                plane_points[:, 0],
                plane_points[:, 1],
                s=mesh_point_size,
                facecolors=mesh_facecolors,
                edgecolors=mesh_edgecolors,
                linewidths=mesh_linewidth,
                zorder=3,
            )
        ax.set_title(spec.title)
        ax.set_xlabel(_axis_label(spec.axis_labels[0], axis_unit))
        ax.set_ylabel(_axis_label(spec.axis_labels[1], axis_unit))
        ax.set_xlim(float(plot_bbox[0, 0]), float(plot_bbox[1, 0]))
        ax.set_ylim(float(plot_bbox[0, 1]), float(plot_bbox[1, 1]))
        ax.set_aspect("equal")

    if created_axes:
        for ax in axes_array[len(plane_specs):]:
            ax.axis("off")

    fig.suptitle(suptitle, fontsize=16)
    fig.legend(
        handles=_legend_handles(display_names),
        loc="outside lower center",
        ncol=min(4, len(display_names)),
        frameon=False,
    )
    return fig, axes_array

plot_scalar_field_on_planes(problem, scalar_values, plane_specs, *, axes=None, ncols=None, figsize=(13.0, 11.0), suptitle='Scalar field cuts', colorbar_label='Value', cmap='coolwarm', symmetric=True, vmin=None, vmax=None, tile_edgecolors='none', tile_linewidth=0.0, show_volume_edges=False, coordinate_scale=1.0, axis_unit=None)

Plot one scalar field on a set of mesh-aligned 2D plane cuts.

Source code in tempura/plotting/meshing.py

def plot_scalar_field_on_planes(
    problem: Any,
    scalar_values: np.ndarray,
    plane_specs: Sequence[PlaneSpec],
    *,
    axes: Sequence[plt.Axes] | None = None,
    ncols: int | None = None,
    figsize: tuple[float, float] = (13.0, 11.0),
    suptitle: str = "Scalar field cuts",
    colorbar_label: str = "Value",
    cmap: str = "coolwarm",
    symmetric: bool = True,
    vmin: float | None = None,
    vmax: float | None = None,
    tile_edgecolors: str = "none",
    tile_linewidth: float = 0.0,
    show_volume_edges: bool = False,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
) -> tuple[plt.Figure, np.ndarray]:
    """Plot one scalar field on a set of mesh-aligned 2D plane cuts."""
    if not plane_specs:
        raise ValueError("plane_specs must not be empty.")
    if ncols is None:
        ncols = 2 if len(plane_specs) > 1 else 1
    if ncols <= 0:
        raise ValueError(f"ncols must be positive; got {ncols}.")

    coordinates = np.asarray(problem.coordinates, dtype=float)
    scalar_values = np.asarray(scalar_values, dtype=float)
    if scalar_values.shape != (coordinates.shape[0],):
        raise ValueError(
            "scalar_values must have shape (npoints,); "
            f"got {scalar_values.shape} for coordinates {coordinates.shape}."
        )

    volume_bounds = _volume_bounds_lookups(coordinates)
    panel_payloads: list[tuple[PlaneSpec, np.ndarray, np.ndarray]] = []
    panel_values: list[np.ndarray] = []
    for spec in plane_specs:
        indices = _exact_plane_site_indices(
            coordinates,
            spec=spec,
        )
        if indices.size == 0:
            indices = _plane_site_indices(
                coordinates,
                spec=spec,
                bounds_lookup=volume_bounds[spec.plane_axis],
            )
        sites = coordinates[indices] if indices.size else np.empty((0, 3), dtype=float)
        values = scalar_values[indices] if indices.size else np.empty((0,), dtype=float)
        panel_payloads.append((spec, sites, values))
        if values.size:
            panel_values.append(values)

    if panel_values:
        flattened = np.concatenate(panel_values)
        if symmetric:
            color_limit = float(np.max(np.abs(flattened)))
            if np.isclose(color_limit, 0.0):
                color_limit = 1.0
            resolved_vmin = -color_limit if vmin is None else float(vmin)
            resolved_vmax = color_limit if vmax is None else float(vmax)
        else:
            resolved_vmin = float(np.min(flattened) if vmin is None else vmin)
            resolved_vmax = float(np.max(flattened) if vmax is None else vmax)
            if np.isclose(resolved_vmin, resolved_vmax):
                resolved_vmin -= 1.0
                resolved_vmax += 1.0
    else:
        resolved_vmin = -1.0 if vmin is None else float(vmin)
        resolved_vmax = 1.0 if vmax is None else float(vmax)

    cmap_obj = plt.cm.get_cmap(cmap)
    norm = Normalize(vmin=resolved_vmin, vmax=resolved_vmax)
    if axes is None:
        nrows = int(np.ceil(len(plane_specs) / ncols))
        fig, axes_obj = plt.subplots(nrows, ncols, figsize=figsize, constrained_layout=True)
        axes_array = np.atleast_1d(axes_obj).ravel()
        created_axes = True
    else:
        axes_array = np.atleast_1d(np.asarray(list(axes), dtype=object)).ravel()
        if len(axes_array) < len(plane_specs):
            raise ValueError(
                f"axes must contain at least {len(plane_specs)} axes; got {len(axes_array)}."
            )
        fig = axes_array[0].figure
        created_axes = False

    for ax, (spec, plane_sites, values) in zip(axes_array, panel_payloads, strict=False):
        plot_bbox = np.asarray(spec.plot_region.bbox, dtype=float) * float(coordinate_scale)
        plane_grid = _scalar_plane_grid(plane_sites, values, spec=spec)
        if plane_grid is not None:
            x_values, y_values, display_values = plane_grid
            ax.contourf(
                x_values * float(coordinate_scale),
                y_values * float(coordinate_scale),
                display_values,
                levels=np.linspace(resolved_vmin, resolved_vmax, 33),
                cmap=cmap_obj,
                norm=norm,
                antialiased=False,
                zorder=1,
            )
        if show_volume_edges:
            volume_edges = _volume_edges(
                plane_sites,
                spec=spec,
                coordinate_scale=coordinate_scale,
                volume_bounds=volume_bounds,
                tile_edgecolors=tile_edgecolors,
                tile_linewidth=tile_linewidth,
            )
            if volume_edges is not None:
                ax.add_collection(volume_edges)
        ax.set_title(spec.title)
        ax.set_xlabel(_axis_label(spec.axis_labels[0], axis_unit))
        ax.set_ylabel(_axis_label(spec.axis_labels[1], axis_unit))
        ax.set_xlim(float(plot_bbox[0, 0]), float(plot_bbox[1, 0]))
        ax.set_ylim(float(plot_bbox[0, 1]), float(plot_bbox[1, 1]))
        ax.set_aspect("equal")

    if created_axes:
        for ax in axes_array[len(plane_specs):]:
            ax.axis("off")

    sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap_obj)
    sm.set_array([])
    fig.suptitle(suptitle, fontsize=16)
    fig.colorbar(
        sm,
        ax=list(axes_array[:len(plane_specs)]),
        label=colorbar_label,
        location="right",
        pad=0.03,
        fraction=0.05,
    )
    return fig, axes_array

save_problem_regions_with_mesh(problem, region_shapes, plane_specs, output_path, *, region_names=None, region_display_names=None, ncols=None, figsize=(13.0, 11.0), suptitle='Pescado region cuts with overlaid mesh nodes', dpi=220, fill_mode='regions', background_alpha=0.88, tile_edgecolors='white', tile_linewidth=0.6, show_volume_edges=False, show_mesh_points=True, mesh_point_size=12.0, mesh_facecolors='white', mesh_edgecolors='black', mesh_linewidth=0.35, coordinate_scale=1.0, axis_unit=None)

Render region slices with overlaid mesh nodes and save the figure.

This is a convenience wrapper around plot_problem_regions_with_mesh(...) that also creates the parent directory and closes the Matplotlib figure after saving.

Source code in tempura/plotting/meshing.py

def save_problem_regions_with_mesh(
    problem: Any,
    region_shapes: RegionShapeMap,
    plane_specs: Sequence[PlaneSpec],
    output_path: str | Path,
    *,
    region_names: Sequence[str] | None = None,
    region_display_names: RegionDisplayMap | None = None,
    ncols: int | None = None,
    figsize: tuple[float, float] = (13.0, 11.0),
    suptitle: str = "Pescado region cuts with overlaid mesh nodes",
    dpi: int = 220,
    fill_mode: str = "regions",
    background_alpha: float = 0.88,
    tile_edgecolors: str = "white",
    tile_linewidth: float = 0.6,
    show_volume_edges: bool = False,
    show_mesh_points: bool = True,
    mesh_point_size: float = 12.0,
    mesh_facecolors: str = "white",
    mesh_edgecolors: str = "black",
    mesh_linewidth: float = 0.35,
    coordinate_scale: float = 1.0,
    axis_unit: str | None = None,
) -> Path:
    """Render region slices with overlaid mesh nodes and save the figure.

    This is a convenience wrapper around ``plot_problem_regions_with_mesh(...)``
    that also creates the parent directory and closes the Matplotlib figure
    after saving.
    """
    output_path = Path(output_path)
    figure, _ = plot_problem_regions_with_mesh(
        problem,
        region_shapes,
        plane_specs,
        region_names=region_names,
        region_display_names=region_display_names,
        ncols=ncols,
        figsize=figsize,
        suptitle=suptitle,
        fill_mode=fill_mode,
        background_alpha=background_alpha,
        tile_edgecolors=tile_edgecolors,
        tile_linewidth=tile_linewidth,
        show_volume_edges=show_volume_edges,
        show_mesh_points=show_mesh_points,
        mesh_point_size=mesh_point_size,
        mesh_facecolors=mesh_facecolors,
        mesh_edgecolors=mesh_edgecolors,
        mesh_linewidth=mesh_linewidth,
        coordinate_scale=coordinate_scale,
        axis_unit=axis_unit,
    )
    output_path.parent.mkdir(parents=True, exist_ok=True)
    figure.savefig(output_path, dpi=dpi)
    plt.close(figure)
    return output_path

transpose_plane_spec(spec, *, title=None)

Return one plane spec with the displayed axes swapped.

Source code in tempura/plotting/meshing.py

def transpose_plane_spec(spec: PlaneSpec, *, title: str | None = None) -> PlaneSpec:
    """Return one plane spec with the displayed axes swapped."""
    bbox = np.asarray(spec.plot_region.bbox, dtype=float)
    lower_left = np.asarray(bbox[0, ::-1], dtype=float)
    size = np.asarray((bbox[1] - bbox[0])[::-1], dtype=float)
    return replace(
        spec,
        title=spec.title if title is None else title,
        plot_region=shapes.Box(lower_left=lower_left, size=size),
        axis_indices=tuple(reversed(spec.axis_indices)),
        axis_labels=tuple(reversed(spec.axis_labels)),
    )

Viewer

Viewer export helpers for serialized device scenes.

DeviceViewExportManifest dataclass

Paths written by :func:export_device_view.

Attributes:

Name	Type	Description
scene_name	str	Human-readable scene name stored in scene.json.
scene_json_path	Path	Path to the scene metadata file.
asset_paths	tuple[Path, ...]	Paths to binary array assets referenced by scene.json.

Source code in tempura/viewer/export.py

@dataclass(frozen=True)
class DeviceViewExportManifest:
    """Paths written by :func:`export_device_view`.

    Attributes:
        scene_name: Human-readable scene name stored in ``scene.json``.
        scene_json_path: Path to the scene metadata file.
        asset_paths: Paths to binary array assets referenced by ``scene.json``.
    """

    scene_name: str
    scene_json_path: Path
    asset_paths: tuple[Path, ...]

export_device_view(device, out_dir, *, scene_name='device', colors=None, overwrite=True, verbose=False)

Export a device into files consumed by the TypeScript 3D viewer.

The exporter writes one scene.json metadata file and one or more binary array files storing layer height fields (plus gate stencils).

Parameters:

Name	Type	Description	Default
device	Device	Device to export.	required
out_dir	str | Path	Output directory for the exported scene bundle.	required
scene_name	str	Label stored in the metadata file.	'device'
colors	dict[str, str] | None	Optional per-layer color overrides keyed by layer name.	None
overwrite	bool	If False, refuse to write into a non-empty directory.	True
verbose	bool	If True, print progress messages.	False

Returns:

Type	Description
DeviceViewExportManifest	Manifest describing the files written.

Raises:

Type	Description
ValueError	If the device is empty or a color override is invalid.
FileExistsError	If overwrite is disabled and the destination exists.
RuntimeError	If the device is not finalized or finalized height arrays are missing.

Source code in tempura/viewer/export.py

def export_device_view(
    device: Device,
    out_dir: str | Path,
    *,
    scene_name: str = "device",
    colors: dict[str, str] | None = None,
    overwrite: bool = True,
    verbose: bool = False,
) -> DeviceViewExportManifest:
    """Export a device into files consumed by the TypeScript 3D viewer.

    The exporter writes one ``scene.json`` metadata file and one or more binary
    array files storing layer height fields (plus gate stencils).

    Args:
        device: Device to export.
        out_dir: Output directory for the exported scene bundle.
        scene_name: Label stored in the metadata file.
        colors: Optional per-layer color overrides keyed by layer name.
        overwrite: If ``False``, refuse to write into a non-empty directory.
        verbose: If ``True``, print progress messages.

    Returns:
        Manifest describing the files written.

    Raises:
        ValueError: If the device is empty or a color override is invalid.
        FileExistsError: If ``overwrite`` is disabled and the destination exists.
        RuntimeError: If the device is not finalized or finalized height arrays are missing.
    """

    def _log(message: str) -> None:
        """Print a progress message when verbose logging is enabled."""
        if verbose:
            print(message)

    if not device.layers:
        raise ValueError("Device must contain at least one layer before export.")

    if not device.finalized:
        raise RuntimeError("Device must be finalized before export_device_view().")

    out_path = Path(out_dir)
    if out_path.exists():
        if not out_path.is_dir():
            raise FileExistsError(f"Output path is not a directory: {out_path}")
        if not overwrite and any(out_path.iterdir()):
            raise FileExistsError(f"Output directory already exists and is not empty: {out_path}")
    else:
        out_path.mkdir(parents=True, exist_ok=True)

    arrays_dir = out_path / "arrays"
    arrays_dir.mkdir(parents=True, exist_ok=True)

    color_map = _default_layer_colors(device)
    if colors is not None:
        unknown_layers = sorted(set(colors) - set(device.layers))
        if unknown_layers:
            raise ValueError(f"Unknown layer names in colors override: {unknown_layers}")
        color_map.update(colors)

    asset_paths: list[Path] = []
    layers_payload = []
    for order, layer in enumerate(device.layers.values()):
        _log(f"[export_device_view] Exporting layer {layer.name!r}...")
        layers_payload.append(
            _export_layer(
                layer=layer,
                order=order,
                arrays_dir=arrays_dir,
                color=color_map[layer.name],
                asset_paths=asset_paths,
            )
        )

    metadata = _serialize_scene_metadata(
        device=device,
        scene_name=scene_name,
        layers_payload=layers_payload,
    )
    scene_json_path = out_path / "scene.json"
    scene_json_path.write_text(json.dumps(metadata, indent=2) + "\n")
    _log(f"[export_device_view] Wrote {scene_json_path}")

    return DeviceViewExportManifest(
        scene_name=scene_name,
        scene_json_path=scene_json_path,
        asset_paths=tuple(asset_paths),
    )