electrostatics._pescado_solve — gate-basis orchestration

tempura/electrostatics/_pescado_solve.py orchestrates the gate-basis solve. It reads the region structure back from a finalized pescado problem, hands the matrix to the linear core in solver, and returns one basis potential per gate. It backs the public solve_gate_potentials.

The gate basis

solve_gate_potentials(...) computes one potential per gate. For gate \(i\) it sets that gate to \(1\,\mathrm V\), every other gate to \(0\,\mathrm V\), and the outer boundary to \(0\,\mathrm V\), then solves for the resulting field \(\phi_i\).

Because the problem is linear, any real gate configuration with voltages \(V_i\) is the superposition

\[ \phi(\mathbf r) = \sum_i V_i\,\phi_i(\mathbf r). \]

The expensive work — geometry assembly and the matrix factorization — happens once; every new voltage configuration is then just a weighted sum of the stored basis fields.

How it is organized

The orchestration is deliberately stateless with respect to Tempura: it reconstructs everything it needs from the finalized problem's metadata.

1. get_region_inds(...) recovers the solver indices for the interior, each gate, and the boundary from the finalized problem's region shapes.
2. _prepare_gate_solve(...) builds the shared symmetric Dirichlet system once (see solver) along with the gate ordering taken from the problem's Dirichlet region order.
3. _solve_gate_blocks(...) factors that system once and solves the gate basis in blocks of rhs_block_size columns, so a large gate set never has to materialize every right-hand side at once.

The gate order is read from problem.regions["dirichlet"] (excluding the Boundary entry), which fixes a stable, reproducible order for the returned basis mapping.

When superposition breaks down

Superposition assumes fixed geometry, fixed materials, and a linear response. If the quantum-region response becomes self-consistent, voltage-dependent, or otherwise nonlinear, the basis fields no longer describe the device on their own — hand the finalized problem to Pescado's self-consistent solver directly.

API

_pescado_solve

Gate-solve orchestration for Tempura's active electrostatics pipeline.

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

solve_gate_potentials_components(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 and return the internal solve outputs.

Parameters:

Name	Type	Description	Default
problem	Any	Finalized Pescado problem.	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	str | float	Which realized quantum region z plane to export when the region spans multiple z coordinates.	'midpoint'
verbose	bool	Whether to emit progress logging.	False

Returns:

Type	Description
PotentialMap	Tuple of (basis_potentials, plane). plane is None unless a
QuantumRegionPlaneData | None	quantum region export was requested.

Source code in tempura/electrostatics/_pescado_solve.py

def solve_gate_potentials_components(
    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: str | float = "midpoint",
    verbose: bool = False,
) -> tuple[PotentialMap, QuantumRegionPlaneData | None]:
    """Solve the gate basis problem and return the internal solve outputs.

    Args:
        problem: Finalized Pescado problem.
        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.

    Returns:
        Tuple of ``(basis_potentials, plane)``. ``plane`` is ``None`` unless a
        quantum region export was requested.
    """
    region_shapes = _region_shapes_from_problem(problem)
    _validate_named_regions_have_sites(problem, region_shapes)
    gate_names = _gate_names_from_problem(problem)
    _log_if_verbose(verbose, f"[solve_gate_potentials] Gates: {gate_names}")
    if not gate_names:
        _log_if_verbose(
            verbose,
            "[solve_gate_potentials] No gates provided; returning empty result.",
        )
        return {}, None

    setup = _prepare_gate_solve(
        problem,
        verbose=verbose,
    )
    gate_names = setup["gate_names"]
    results = _solve_gate_blocks(
        setup,
        gate_names,
        charge=charge,
        dtype=dtype,
        blr=blr,
        eps_blr=eps_blr,
        rhs_block_size=rhs_block_size,
        verbose=verbose,
    )

    plane = None
    if save_quantum_region_potential is not None:
        plane = _export_quantum_region_potential_bundle(
            results,
            problem,
            gate_names,
            save_quantum_region_potential=save_quantum_region_potential,
            quantum_region_name=quantum_region_name,
            plane_selection=quantum_region_plane_selection,
            verbose=verbose,
        )

    _log_if_verbose(verbose, "[solve_gate_potentials] Done.")
    return results, plane