electrostatics.pescado_wrapper — the public solve facade

tempura/electrostatics/pescado_wrapper.py is the thin, stable facade over the electrostatics pipeline. It exposes the two functions most users call — build_problem(...) and solve_gate_potentials(...) — and re-exports the quantum-region extraction helpers. The heavy lifting lives in the internal modules it delegates to; this file's job is to keep a small, documented surface.

The two-step solve contract

The workflow is intentionally explicit about where geometry ends and the solve begins:

problem_builder = build_problem(device)        # assemble the meshed problem
problem = problem_builder.finalized()           # freeze it into a pescado problem
basis_potentials = solve_gate_potentials(problem)  # one basis field per gate

build_problem(...) returns a ProblemBuilder (see _pescado_build); you finalize it yourself, then pass the finalized problem to solve_gate_potentials(...) (see _pescado_solve). The facade deliberately does not finalize for you, so the finalize step — and the chance to inspect or hand the problem elsewhere — stays visible. Passing an un-finalized builder to the solve raises a clear error.

What you get back

solve_gate_potentials(...) returns a mapping from gate name to the solved basis-potential vector. Because the solve is linear, any voltage configuration is the superposition \(\phi = \sum_i V_i\,\phi_i\) of those basis fields — the factorization is reused, not recomputed, for each new voltage setting.

Handing off to a self-consistent solver

For a nonlinear or voltage-dependent quantum-region response, the basis fields no longer describe the device on their own. The finalized problem is itself the handoff object — it carries the coordinates, region shapes, sparse-vector helpers, and system matrix a self-consistent solver needs:

from pescado.self_consistent import solver as tf_backend

Tempura does not wrap that solver as part of its core API; it hands you the finalized problem and steps out of the way.

API

pescado_wrapper

Compatibility facade for Tempura's active electrostatics pipeline.

QuantumRegionPlaneData dataclass

Rectangular quantum region plane data extracted from a solved problem.

Attributes:

Name	Type	Description
sorted_indices	NDArray[int_]	Solver-vector indices sorted lexicographically by (y, x) on the selected plane.
coords	NDArray[float64]	Sorted (x, y, z) coordinates on the selected quantum region plane.
x_values	NDArray[float64]	Unique x coordinates spanning the rectangular plane.
y_values	NDArray[float64]	Unique y coordinates spanning the rectangular plane.
z_value	float	Shared z coordinate of the selected plane.
nx	int	Number of x samples in the plane.
ny	int	Number of y samples in the plane.
region_mask	NDArray[bool_]	(ny, nx) boolean mask marking which plane cells are actually owned by the quantum region. For a full-footprint region every cell is True; for a patterned region the plane is a complete rectangle (so basis vectors reshape cleanly) but only the stenciled cells are True, so callers never treat the bounding-box filler cells as active material.

Source code in tempura/electrostatics/_quantum_region.py

@dataclass(frozen=True)
class QuantumRegionPlaneData:
    """Rectangular quantum region plane data extracted from a solved problem.

    Attributes:
        sorted_indices: Solver-vector indices sorted lexicographically by
            ``(y, x)`` on the selected plane.
        coords: Sorted ``(x, y, z)`` coordinates on the selected quantum region plane.
        x_values: Unique x coordinates spanning the rectangular plane.
        y_values: Unique y coordinates spanning the rectangular plane.
        z_value: Shared z coordinate of the selected plane.
        nx: Number of x samples in the plane.
        ny: Number of y samples in the plane.
        region_mask: ``(ny, nx)`` boolean mask marking which plane cells are
            actually owned by the quantum region. For a full-footprint region
            every cell is ``True``; for a **patterned** region the plane is a
            complete rectangle (so basis vectors reshape cleanly) but only the
            stenciled cells are ``True``, so callers never treat the
            bounding-box filler cells as active material.
    """

    sorted_indices: NDArray[np.int_]
    coords: NDArray[np.float64]
    x_values: NDArray[np.float64]
    y_values: NDArray[np.float64]
    z_value: float
    nx: int
    ny: int
    region_mask: NDArray[np.bool_]

build_problem(device, vacuum_scale=8.0, *, vacuum_resolution_scale=8.0, dirichlet_boundary=True, 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
vacuum_resolution_scale	float	Coarsening factor applied to the automatic outer-vacuum lattice spacing.	8.0
dirichlet_boundary	bool	Whether to add a thin outer shell of zero-volt Dirichlet sites around the simulation box. When False, the outer shell is omitted.	True
verbose	bool	Whether to emit progress logging during assembly.	False

Returns:

Type	Description
ProblemBuilder	Initialized Pescado :class:ProblemBuilder ready to be finalized
ProblemBuilder	explicitly before calling :func:solve_gate_potentials.

Source code in tempura/electrostatics/pescado_wrapper.py

def build_problem(
    device: Device,
    vacuum_scale: float = 8.0,
    *,
    vacuum_resolution_scale: float = 8.0,
    dirichlet_boundary: bool = True,
    verbose: bool = False,
) -> ProblemBuilder:
    """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.
        vacuum_resolution_scale: Coarsening factor applied to the automatic
            outer-vacuum lattice spacing.
        dirichlet_boundary: Whether to add a thin outer shell of zero-volt
            Dirichlet sites around the simulation box. When ``False``, the
            outer shell is omitted.
        verbose: Whether to emit progress logging during assembly.

    Returns:
        Initialized Pescado :class:`ProblemBuilder` ready to be finalized
        explicitly before calling :func:`solve_gate_potentials`.
    """
    problem_builder, _region_shapes, _gate_names, _ = build_problem_components(
        device,
        vacuum_scale=vacuum_scale,
        vacuum_resolution_scale=vacuum_resolution_scale,
        dirichlet_boundary=dirichlet_boundary,
        verbose=verbose,
    )
    return problem_builder

extract_quantum_region_basis(basis_potentials, gate_names, plane)

Reshape solved gate basis vectors onto a rectangular quantum region 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	QuantumRegionPlaneData	Extracted quantum region plane metadata and sorted solver indices.	required

Returns:

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

Source code in tempura/electrostatics/_quantum_region.py

def extract_quantum_region_basis(
    basis_potentials: Mapping[str, np.ndarray],
    gate_names: list[str],
    plane: QuantumRegionPlaneData,
) -> dict[str, np.ndarray]:
    """Reshape solved gate basis vectors onto a rectangular quantum region plane.

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

    Returns:
        Mapping of gate name to ``(ny, nx)`` quantum region 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
    }

extract_quantum_region_plane(problem, region_shapes, *, region_name='quantum_region', plane_selection='midpoint')

Return one rectangular quantum region 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.	required
region_name	str	Region to extract from region_shapes.	'quantum_region'
plane_selection	PlaneSelection	Which z plane to export when the realized quantum region spans multiple z coordinates. Use "lower", "midpoint", "upper", or an explicit z value.	'midpoint'

Returns:

Type	Description
QuantumRegionPlaneData	Sorted rectangular plane data for the selected z slice across the full
QuantumRegionPlaneData	device XY span.

Raises:

Type	Description
ValueError	If the quantum region is missing, empty, or the selected sampling slab does not form a rectangular plane.

Source code in tempura/electrostatics/_quantum_region.py

def extract_quantum_region_plane(
    problem: Any,
    region_shapes: Mapping[str, Any],
    *,
    region_name: str = "quantum_region",
    plane_selection: PlaneSelection = "midpoint",
) -> QuantumRegionPlaneData:
    """Return one rectangular quantum region plane extracted from ``problem``.

    Args:
        problem: Finalized Pescado problem exposing ``coordinates`` and
            ``points_inside``.
        region_shapes: Mapping of region names to shapes.
        region_name: Region to extract from ``region_shapes``.
        plane_selection: Which z plane to export when the realized quantum region spans
            multiple z coordinates. Use ``"lower"``, ``"midpoint"``,
            ``"upper"``, or an explicit z value.

    Returns:
        Sorted rectangular plane data for the selected z slice across the full
        device XY span.

    Raises:
        ValueError: If the quantum region is missing, empty, or the selected
            sampling slab does not form a rectangular plane.
    """
    try:
        region_shape = region_shapes[region_name]
    except KeyError as exc:
        raise ValueError(
            "quantum region export requires a matching region_shapes entry; "
            f"missing {region_name!r}."
        ) from exc

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

    z_vals = np.unique(np.round(coords_quantum_region[:, 2], 12))
    z_target = _resolve_plane_z_value(z_vals, plane_selection=plane_selection)
    # Sample the quantum region's own XY/z extent to recover a rectangular
    # plane before slicing to the requested z level. Using the region's own
    # bounding box (rather than a symmetric max over every region) keeps the
    # slab off neighboring coarser layers and works for devices that are not
    # centered on the origin.
    slab = _sampling_slab_for_quantum_region_plane(region_shape)
    indices_quantum_region = np.asarray(problem.points_inside(slab), dtype=int)
    coords_quantum_region = np.asarray(
        problem.coordinates[indices_quantum_region],
        dtype=float,
    )
    if coords_quantum_region.size == 0:
        raise ValueError(
            "No solver coordinates were found inside the quantum region sampling slab."
        )

    plane_mask = np.isclose(coords_quantum_region[:, 2], z_target, atol=1e-12, rtol=0.0)
    if not np.any(plane_mask):
        raise ValueError(
            "The requested quantum region plane does not match any realized "
            f"z coordinates; requested={z_target}, available={z_vals.tolist()}."
        )
    coords_quantum_region = coords_quantum_region[plane_mask]
    indices_quantum_region = indices_quantum_region[plane_mask]

    # Sort by y first, then x, so reshape(ny, nx) produces row-major arrays
    # matching plotting/image conventions.
    order = np.lexsort((coords_quantum_region[:, 0], coords_quantum_region[:, 1]))
    coords_quantum_region = coords_quantum_region[order]
    sorted_indices_quantum_region = indices_quantum_region[order]

    # Rectangular-grid validation catches malformed or partially meshed active
    # regions before callers reshape solver vectors into plane arrays.
    x_values = np.unique(np.round(coords_quantum_region[:, 0], 12))
    y_values = np.unique(np.round(coords_quantum_region[:, 1], 12))
    nx = len(x_values)
    ny = len(y_values)
    if coords_quantum_region.shape[0] != nx * ny:
        raise ValueError(
            "The extracted quantum region coordinates do not form a rectangular grid."
        )

    x_grid = coords_quantum_region[:, 0].reshape(ny, nx)
    y_grid = coords_quantum_region[:, 1].reshape(ny, nx)
    if not np.allclose(x_grid, np.broadcast_to(x_values, (ny, nx))):
        raise ValueError(
            "The extracted quantum region 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 quantum region y coordinates do not form a rectangular grid."
        )

    # The plane is the region's full bounding box, so a patterned region leaves
    # non-owned filler cells inside it. Mark which cells the region actually
    # owns (the same membership test as ``points_inside``) so callers can mask
    # out the filler instead of treating it as active material.
    region_mask = np.asarray(region_shape(coords_quantum_region), dtype=bool).reshape(
        ny, nx
    )

    return QuantumRegionPlaneData(
        sorted_indices=sorted_indices_quantum_region,
        coords=coords_quantum_region,
        x_values=np.asarray(x_values, dtype=float),
        y_values=np.asarray(y_values, dtype=float),
        z_value=float(coords_quantum_region[0, 2]),
        nx=nx,
        ny=ny,
        region_mask=region_mask,
    )

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

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

Parameters:

Name	Type	Description	Default
pescado_problem_finalized	Any	Finalized Pescado problem exposing region index arrays plus points_inside.	required
gate_shapes	RegionMap	Mapping of gate name to solver-owned gate region shape.	required
boundary_shape	Shape | None	Solver-owned outer Dirichlet boundary shape, or None when no outer boundary is present.	required
boundary_name	str	Key used for the boundary entry in the returned mapping.	'Boundary'

Returns:

Type	Description
IndexMap	Mapping containing "Interior" and one entry per gate in
IndexMap	gate_shapes. When boundary_shape is not None, the mapping
IndexMap	also contains one entry for boundary_name.

Source code in tempura/electrostatics/_pescado_solve.py

def get_region_inds(
    pescado_problem_finalized: Any,
    gate_shapes: RegionMap,
    boundary_shape: shapes.Shape | None,
    *,
    boundary_name: str = "Boundary",
) -> IndexMap:
    """Return solver indices for interior, gates, and the outer boundary.

    Args:
        pescado_problem_finalized: Finalized Pescado problem exposing region
            index arrays plus ``points_inside``.
        gate_shapes: Mapping of gate name to solver-owned gate region shape.
        boundary_shape: Solver-owned outer Dirichlet boundary shape, or
            ``None`` when no outer boundary is present.
        boundary_name: Key used for the boundary entry in the returned mapping.

    Returns:
        Mapping containing ``"Interior"`` and one entry per gate in
        ``gate_shapes``. When ``boundary_shape`` is not ``None``, the mapping
        also contains one entry for ``boundary_name``.
    """
    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)

    if boundary_shape is not None:
        region_inds[boundary_name] = pescado_problem_finalized.points_inside(
            boundary_shape
        )
    return region_inds

project_gate_footprints_to_plane(gate_shapes, plane)

Project gate XY footprints onto the exported rectangular quantum region plane.

Each gate is sampled at the midpoint of its z extent so the resulting mask matches the gate's in-plane footprint while remaining aligned to the quantum region export grid.

Parameters:

Name	Type	Description	Default
gate_shapes	Mapping[str, Any]	Mapping of gate name to callable solver-owned gate shape.	required
plane	QuantumRegionPlaneData	Rectangular quantum region plane onto which the footprints should be sampled.	required

Returns:

Type	Description
dict[str, ndarray]	Mapping of gate name to (ny, nx) boolean footprint mask aligned to
dict[str, ndarray]	plane.

Source code in tempura/electrostatics/_quantum_region.py

def project_gate_footprints_to_plane(
    gate_shapes: Mapping[str, Any],
    plane: QuantumRegionPlaneData,
) -> dict[str, np.ndarray]:
    """Project gate XY footprints onto the exported rectangular quantum region plane.

    Each gate is sampled at the midpoint of its z extent so the resulting mask
    matches the gate's in-plane footprint while remaining aligned to the quantum region
    export grid.

    Args:
        gate_shapes: Mapping of gate name to callable solver-owned gate shape.
        plane: Rectangular quantum region plane onto which the footprints should be
            sampled.

    Returns:
        Mapping of gate name to ``(ny, nx)`` boolean footprint mask aligned to
        ``plane``.
    """
    if not gate_shapes:
        return {}

    coords = np.asarray(plane.coords, dtype=float)
    footprints: dict[str, np.ndarray] = {}
    for gate_name, gate_shape in gate_shapes.items():
        # The footprint projection only needs the gate's XY shape. Sampling at
        # the gate midpoint in z avoids coupling the mask to the quantum-region
        # plane height.
        bbox = getattr(gate_shape, "bbox", None)
        if bbox is None or not callable(gate_shape):
            continue

        bbox_arr = np.asarray(bbox, dtype=float)
        if bbox_arr.shape != (2, 3):
            continue
        z_sample = float(0.5 * (bbox_arr[0, 2] + bbox_arr[1, 2]))
        sample_points = coords.copy()
        sample_points[:, 2] = z_sample

        mask = np.asarray(gate_shape(sample_points), dtype=bool).reshape(-1)
        if mask.shape[0] != coords.shape[0]:
            continue
        footprints[gate_name] = mask.reshape(plane.ny, plane.nx)

    return footprints

solve_gate_potentials(problem, *, charge=None, dtype=DTYPE, blr=BLR, eps_blr=EPS_BLR, rhs_block_size=8, save_quantum_region_potential=None, quantum_region_name='quantum_region', quantum_region_plane_selection='midpoint', verbose=False)

Solve the gate basis problem with one reused factorization.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem returned by build_problem(...).finalized().	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_quantum_region_potential	str | Path | None	Optional output directory for static quantum region export.	None
quantum_region_name	str	Region name to export when saving a quantum region basis bundle.	'quantum_region'
quantum_region_plane_selection	PlaneSelection	Which realized quantum region 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	Mapping from gate name to solved basis-potential vector.

Source code in tempura/electrostatics/pescado_wrapper.py

def solve_gate_potentials(
    problem: Any,
    *,
    charge: np.ndarray | None = None,
    dtype: DTypeLike = DTYPE,
    blr: bool = BLR,
    eps_blr: float = EPS_BLR,
    rhs_block_size: int = 8,
    save_quantum_region_potential: str | Path | None = None,
    quantum_region_name: str = "quantum_region",
    quantum_region_plane_selection: PlaneSelection = "midpoint",
    verbose: bool = False,
) -> PotentialMap:
    """Solve the gate basis problem with one reused factorization.

    Args:
        problem: Finalized Pescado problem returned by
            ``build_problem(...).finalized()``.
        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_quantum_region_potential: Optional output directory for static
            quantum region export.
        quantum_region_name: Region name to export when saving a quantum region
            basis bundle.
        quantum_region_plane_selection: Which realized quantum region z plane to
            export when the region spans multiple z coordinates.
        verbose: Whether to emit progress logging during solving.

    Returns:
        Mapping from gate name to solved basis-potential vector.
    """
    if isinstance(problem, ProblemBuilder):
        raise TypeError(
            "solve_gate_potentials() expects a finalized Pescado problem; call "
            "build_problem(...).finalized() first."
        )
    _region_shapes_from_problem(problem)

    basis_potentials, _plane = solve_gate_potentials_components(
        problem,
        charge=charge,
        dtype=dtype,
        blr=blr,
        eps_blr=eps_blr,
        rhs_block_size=rhs_block_size,
        save_quantum_region_potential=save_quantum_region_potential,
        quantum_region_name=quantum_region_name,
        quantum_region_plane_selection=quantum_region_plane_selection,
        verbose=verbose,
    )
    return basis_potentials