API

Previous

Device model

tempura.electrostatics.device

Device geometry helpers for building meshing regions.

Device

Container for device layers and derived height/shape fields.

Attributes:

Name Type Description
L

Device length.
W

Device width.
dx

Grid spacing in x.
dy

Grid spacing in y.
min_grid_resolution

Global minimum mesh resolution used by electrostatics.
grid_shape

Array shape (rows, cols) for masks.
layers dict[str, Layer]

Layers in stacking order.
Source code in tempura/electrostatics/device.py 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
class Device:
    """Container for device layers and derived height/shape fields.

    Attributes:
        L: Device length.
        W: Device width.
        dx: Grid spacing in x.
        dy: Grid spacing in y.
        min_grid_resolution: Global minimum mesh resolution used by electrostatics.
        grid_shape: Array shape (rows, cols) for masks.
        layers: Layers in stacking order.
    """

    MIN_GRID_RESOLUTION = 1.0

    def __init__(self, device_size: tuple[float, float, float, float]):
        """Create a device with (L, W, dx, dy).

        Args:
            device_size: Tuple of (L, W, dx, dy).
        """
        self.L, self.W, self.dx, self.dy = device_size
        self.min_grid_resolution = float(self.MIN_GRID_RESOLUTION)
        self.grid_shape = (
            int(round(self.W / self.dy)),
            int(round(self.L / self.dx)),
        )
        self.layers: dict[str, Layer] = {}
        self._finalized = False

    def add_layer(self, layer: Layer) -> None:
        """Add a layer (order matters for stacking).

        Args:
            layer: Layer instance to add.
        """
        self.layers[layer.name] = layer
        self._finalized = False

    def finalize(self) -> None:
        """Compute height fields and shapes for all layers."""
        self._build_height_fields()
        self._build_shapes()
        self._finalized = True

    def _iter_deposition_batches(self) -> Iterator[list[Layer]]:
        """Yield layers grouped by deposition step.

        Non-gate layers are emitted one at a time. Consecutive gates sharing
        the same deposition group are emitted together so they can occupy the
        same physical z-level during finalization.
        """
        layers = list(self.layers.values())
        index = 0

        while index < len(layers):
            layer = layers[index]
            if not isinstance(layer, Gate):
                yield [layer]
                index += 1
                continue

            gate_batch: list[Layer] = [layer]
            current_gate = layer
            index += 1
            while index < len(layers):
                next_layer = layers[index]
                if not isinstance(next_layer, Gate):
                    break
                if not self._same_gate_deposition_group(current_gate, next_layer):
                    break
                gate_batch.append(next_layer)
                current_gate = next_layer
                index += 1
            yield gate_batch

    def _build_height_fields(self) -> None:
        """Compute per-layer height fields from masks and thicknesses."""
        batches = list(self._iter_deposition_batches())
        h = np.zeros(self.grid_shape, dtype=float)
        fields_by_id: dict[int, tuple[np.ndarray, np.ndarray]] = {}

        for batch_index, batch in enumerate(batches):
            layer = batch[0]
            next_layer = (
                batches[batch_index + 1][0] if batch_index + 1 < len(batches) else None
            )
            if isinstance(layer, Gate):
                gate_batch = [gate for gate in batch if isinstance(gate, Gate)]
                h_before = h.copy()
                group_masks = [self._layer_mask(gate) for gate in gate_batch]

                for gate, mask in zip(gate_batch, group_masks, strict=False):
                    gate_thickness = float(gate.thickness)
                    gate_base = overhang_base(
                        h_before, thickness=gate_thickness, dx=self.dx
                    )
                    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] = gate_base[mask] + gate_thickness
                    fields_by_id[id(gate)] = (z_bot, z_top)

                # Gate batches only change the running topography once a
                # following dielectric buries the metal.
                if next_layer is not None and isinstance(next_layer, Dielectric):
                    group_union_mask = np.logical_or.reduce(
                        np.stack(group_masks, axis=0)
                    )
                    min_thickness = min(float(gate.thickness) for gate in gate_batch)
                    group_base = overhang_base(
                        h_before, thickness=min_thickness, dx=self.dx
                    )
                    z_top_group = np.full_like(h_before, np.nan, dtype=float)
                    z_top_group[group_union_mask] = (
                        group_base[group_union_mask] + min_thickness
                    )
                    h = np.maximum(h_before, np.nan_to_num(z_top_group, nan=-np.inf))
                else:
                    h = h_before
                continue

            mask = self._layer_mask(layer)
            thickness = float(layer.thickness)
            base = overhang_base(h, thickness=thickness, dx=self.dx)
            z_bot = np.full_like(h, np.nan, dtype=float)
            z_top = np.full_like(h, np.nan, dtype=float)
            z_bot[mask] = h[mask]
            z_top[mask] = base[mask] + thickness
            h = np.maximum(h, np.nan_to_num(z_top, nan=-np.inf))
            fields_by_id[id(layer)] = (z_bot, z_top)

        for layer in self.layers.values():
            z_bot, z_top = fields_by_id[id(layer)]
            # Keep the sampled arrays so downstream exporters can write the
            # exact finalized geometry without re-sampling the callables.
            layer.z_bot_arr = z_bot
            layer.z_top_arr = z_top
            layer.z_bot_fn = make_height_fn(z_bot, self.L, self.W)
            layer.z_top_fn = make_height_fn(z_top, self.L, self.W)

    def _build_shapes(self) -> None:
        """Build geometry shapes for each layer."""
        z0 = 0.0
        batches = list(self._iter_deposition_batches())

        for batch_index, batch in enumerate(batches):
            layer = batch[0]
            next_layer = (
                batches[batch_index + 1][0] if batch_index + 1 < len(batches) else None
            )
            if isinstance(layer, Gate):
                gate_batch = [gate for gate in batch if isinstance(gate, Gate)]

                # Consecutive gates in the same deposition batch share z0.
                for gate in gate_batch:
                    gate._build_shape_from_height(self.L, self.W, z0)

                if next_layer is not None and isinstance(next_layer, Dielectric):
                    # Deposit dielectric above the thinnest gate of the batch.
                    z0 += min(float(gate.thickness) for gate in gate_batch)
                continue

            # Non-gate layers keep standard stacking behavior.
            layer._build_shape_from_height(self.L, self.W, z0)
            z0 += layer.thickness

    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.astype(bool)
        # Dielectrics and 2DEG layers always cover the full XY domain.
        return np.ones(self.grid_shape, dtype=bool)

    @staticmethod
    def _same_gate_deposition_group(lhs: Gate, rhs: Gate) -> bool:
        """Return whether two gates belong to the same deposition batch."""
        if lhs.deposition_group is None and rhs.deposition_group is None:
            return True
        if lhs.deposition_group is None or rhs.deposition_group is None:
            return False
        return lhs.deposition_group == rhs.deposition_group

__init__(device_size)

Create a device with (L, W, dx, dy).

Parameters:

Name Type Description Default
device_size tuple[float, float, float, float]

Tuple of (L, W, dx, dy).

 required
Source code in tempura/electrostatics/device.py 
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def __init__(self, device_size: tuple[float, float, float, float]):
    """Create a device with (L, W, dx, dy).

    Args:
        device_size: Tuple of (L, W, dx, dy).
    """
    self.L, self.W, self.dx, self.dy = device_size
    self.min_grid_resolution = float(self.MIN_GRID_RESOLUTION)
    self.grid_shape = (
        int(round(self.W / self.dy)),
        int(round(self.L / self.dx)),
    )
    self.layers: dict[str, Layer] = {}
    self._finalized = False

add_layer(layer)

Add a layer (order matters for stacking).

Parameters:

Name Type Description Default
layer Layer

Layer instance to add.

 required
Source code in tempura/electrostatics/device.py 
178
179
180
181
182
183
184
185
def add_layer(self, layer: Layer) -> None:
    """Add a layer (order matters for stacking).

    Args:
        layer: Layer instance to add.
    """
    self.layers[layer.name] = layer
    self._finalized = False

finalize()

Compute height fields and shapes for all layers.

Source code in tempura/electrostatics/device.py 
187
188
189
190
191
def finalize(self) -> None:
    """Compute height fields and shapes for all layers."""
    self._build_height_fields()
    self._build_shapes()
    self._finalized = True

Dielectric dataclass

Bases: LayerBase

Dielectric layer covering the full XY domain.

Source code in tempura/electrostatics/device.py 
108
109
110
111
112
@dataclass
class Dielectric(LayerBase):
    """Dielectric layer covering the full XY domain."""

    permitivity: float

Gate dataclass

Bases: LayerBase

Metallic gate layer.

Requires a stencil and does not define a dielectric constant.

Source code in tempura/electrostatics/device.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
@dataclass
class Gate(LayerBase):
    """Metallic gate layer.

    Requires a stencil and does not define a dielectric constant.
    """

    stencil: np.ndarray

    def __post_init__(self) -> None:
        """Validate required attributes after dataclass init."""
        if self.stencil is None:
            raise ValueError("Gate requires a stencil.")

__post_init__()

Validate required attributes after dataclass init.

Source code in tempura/electrostatics/device.py 
102
103
104
105
def __post_init__(self) -> None:
    """Validate required attributes after dataclass init."""
    if self.stencil is None:
        raise ValueError("Gate requires a stencil.")

LayerBase dataclass

Shared data and helpers for device layers.

Attributes:

Name Type Description
name str

Layer name.
resolution ndarray

Mesh resolution for this layer.
thickness float

Layer thickness.
deposition_group str | None

Optional deposition group identifier used to keep multiple gates at the same physical z-level.
z_bot_arr ndarray | None

Cached bottom height field on the device grid.
z_top_arr ndarray | None

Cached top height field on the device grid.
z_bot_fn Callable[[float, float], float] | None

Function returning bottom height at (x, y).
z_top_fn Callable[[float, float], float] | None

Function returning top height at (x, y).
shape Shape | None

Cached geometry shape.
Source code in tempura/electrostatics/device.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
@dataclass
class LayerBase:
    """Shared data and helpers for device layers.

    Attributes:
        name: Layer name.
        resolution: Mesh resolution for this layer.
        thickness: Layer thickness.
        deposition_group: Optional deposition group identifier used to keep
            multiple gates at the same physical z-level.
        z_bot_arr: Cached bottom height field on the device grid.
        z_top_arr: Cached top height field on the device grid.
        z_bot_fn: Function returning bottom height at (x, y).
        z_top_fn: Function returning top height at (x, y).
        shape: Cached geometry shape.
    """

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

    z_bot_arr: np.ndarray | None = field(init=False, default=None)
    z_top_arr: np.ndarray | None = field(init=False, default=None)
    z_bot_fn: Callable[[float, float], float] | None = field(init=False, default=None)
    z_top_fn: Callable[[float, float], float] | None = field(init=False, default=None)
    shape: shapes.Shape | None = field(init=False, default=None)

    def box(self, L: float, W: float, z0: float) -> shapes.Box:
        """Return a simple box covering the full device footprint.

        Args:
            L: Device length.
            W: Device width.
            z0: Bottom z coordinate.

        Returns:
            A box shape covering the device footprint.
        """
        return shapes.Box(
            lower_left=np.array([-L / 2, -W / 2, z0]),
            size=np.array([L, W, self.thickness]),
        )

    def _build_shape_from_height(self, L: float, W: float, z0: float) -> shapes.Shape:
        """Build a shape constrained by the height fields.

        Args:
            L: Device length.
            W: Device width.
            z0: Bottom z coordinate.

        Returns:
            Geometry shape matching the height fields.

        Raises:
            ValueError: If height functions are not set.
        """
        if self.z_bot_fn is None or self.z_top_fn is None:
            raise ValueError("z_bot_fn and z_top_fn must be set before build_shape().")

        bbox = np.array(
            [[-L / 2, -W / 2, z0], [L / 2, W / 2, z0 + self.thickness]],
            dtype=float,
        )

        def inside_fn(x: np.ndarray) -> np.ndarray:
            # Point-in-volume test using per-(x, y) height bounds.
            xy = x[:, :2]
            z = x[:, 2]
            zb = np.array([self.z_bot_fn(xi, yi) for xi, yi in xy])
            zt = np.array([self.z_top_fn(xi, yi) for xi, yi in xy])
            return (z >= zb) & (z <= zt)

        self.shape = shapes.General(func=inside_fn, bbox=bbox)
        return self.shape

box(L, W, z0)

Return a simple box covering the full device footprint.

Parameters:

Name Type Description Default
L float

Device length.

 required
W float

Device width.

 required
z0 float

Bottom z coordinate.

 required

Returns:

Type Description
Box

A box shape covering the device footprint.
Source code in tempura/electrostatics/device.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def box(self, L: float, W: float, z0: float) -> shapes.Box:
    """Return a simple box covering the full device footprint.

    Args:
        L: Device length.
        W: Device width.
        z0: Bottom z coordinate.

    Returns:
        A box shape covering the device footprint.
    """
    return shapes.Box(
        lower_left=np.array([-L / 2, -W / 2, z0]),
        size=np.array([L, W, self.thickness]),
    )

TwoDEG dataclass

Bases: LayerBase

Two-dimensional electron gas (2DEG) layer.

Attributes:

Name Type Description
permitivity float

Relative permittivity for the region.
boundary_condition Literal['neumann', 'helmholtz', 'flexible']

Boundary condition type for the region.
Source code in tempura/electrostatics/device.py 
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
@dataclass
class TwoDEG(LayerBase):
    """Two-dimensional electron gas (2DEG) layer.

    Attributes:
        permitivity: Relative permittivity for the region.
        boundary_condition: Boundary condition type for the region.
    """

    permitivity: float
    boundary_condition: Literal["neumann", "helmholtz", "flexible"] = "helmholtz"

    def __post_init__(self) -> None:
        """Validate boundary condition choice."""
        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}"
            )

    def _build_shape_from_height(self, L: float, W: float, z0: float) -> shapes.Box:
        """Build the 2DEG region as a simple box spanning the device footprint."""
        self.shape = shapes.Box(
            lower_left=np.array([-L / 2, -W / 2, z0]),
            size=np.array([L, W, self.thickness]),
        )
        return self.shape

__post_init__()

Validate boundary condition choice.

Source code in tempura/electrostatics/device.py 
127
128
129
130
131
132
133
134
def __post_init__(self) -> None:
    """Validate boundary condition choice."""
    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}"
        )

Linear solver

tempura.electrostatics.solver

Sparse linear-system helpers for Tempura electrostatics solves.

This module is intentionally independent from Pescado's geometry assembly. It only deals with constrained sparse linear systems, RHS construction, and MUMPS factorization reuse.

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 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
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

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

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

Source code in tempura/electrostatics/solver.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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."""
    solver = _require_mumps().Context()
    _factorize_solver(
        solver,
        matrix,
        dtype=dtype,
        blr=blr,
        eps_blr=eps_blr,
        symmetric=symmetric,
    )
    return solver

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.

Source code in tempura/electrostatics/solver.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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."""
    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

Pescado wrapper

tempura.electrostatics.pescado_wrapper

Pescado-specific problem assembly for Tempura devices.

build_problem(device, vacuum_scale=8.0, vacuum_resolution_factor=1 / 4, *, verbose=False)

Build a Poisson problem for a layered device.

Creates a vacuum simulation region around the device, refines the mesh for each layer, assigns material properties, and sets boundary conditions.

Parameters:

Name Type Description Default
device Any

Device instance with layers added.

 required
vacuum_scale float

Multiplier for the vacuum region size relative to device.

 8.0
vacuum_resolution_factor float

Factor to scale the vacuum resolution relative to device size.

 1 / 4
verbose bool

If True, print progress messages while assembling the problem.

 False

Returns:

Type Description
tuple[ProblemBuilder, RegionMap, list[str]]

Tuple of (problem_builder, region_shapes, gate_names).
Source code in tempura/electrostatics/pescado_wrapper.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
def build_problem(
    device: Any,
    vacuum_scale: float = 8.0,
    vacuum_resolution_factor: float = 1 / 4,
    *,
    verbose: bool = False,
) -> tuple[ProblemBuilder, RegionMap, list[str]]:
    """Build a Poisson problem for a layered device.

    Creates a vacuum simulation region around the device, refines the mesh
    for each layer, assigns material properties, and sets boundary conditions.

    Args:
        device: Device instance with layers added.
        vacuum_scale: Multiplier for the vacuum region size relative to device.
        vacuum_resolution_factor: Factor to scale the vacuum resolution relative
            to device size.
        verbose: If ``True``, print progress messages while assembling the problem.

    Returns:
        Tuple of (problem_builder, region_shapes, gate_names).
    """
    def _log(msg: str) -> None:
        if verbose:
            print(msg)

    if not device._finalized:
        _log("[build_problem] Finalizing device...")
        device.finalize()
    else:
        _log("[build_problem] Device already finalized.")

    layers = list(device.layers.values())
    if not layers:
        raise ValueError("Device has no layers.")
    for layer in layers:
        if getattr(layer, "shape", None) is None:
            raise RuntimeError(
                f"Layer {layer.name!r} is missing shape after device.finalize()."
            )

    # Initialize the problem builder and compute overall device height.
    pb = ProblemBuilder()
    total_height = max(
        float(np.asarray(layer.shape.bbox, dtype=float)[1, 2]) for layer in layers
    )
    _log(f"[build_problem] Layers: {len(layers)} | z_top_max={total_height}")

    # Define a large vacuum box around the device to serve as the simulation region.
    corner = np.array(
        [
            -vacuum_scale * device.L / 2,
            -vacuum_scale * device.W / 2,
            -vacuum_scale * device.L / 2,
        ]
    )
    size = np.array(
        [
            vacuum_scale * device.L,
            vacuum_scale * device.W,
            vacuum_scale * device.L + total_height,
        ]
    )
    vacuum = shapes.Box(lower_left=corner, size=size)
    _log(f"[build_problem] Vacuum region corner={corner}, size={size}")

    # Use a coarse mesh for the initial vacuum region.
    vacuum_res = device.L * vacuum_resolution_factor
    _log(f"[build_problem] Vacuum mesh resolution={vacuum_res}")
    pattern = patterns.Rectangular.constant(
        element_size=vacuum_res * np.ones(3), center=corner + size / 2
    )
    pb.initialize_mesh(simulation_region=vacuum, pattern=pattern)
    _log("[build_problem] Initialized mesh.")

    # Refine each layer region using the layer-defined resolution.
    bbox_eps = 1e-6
    for layer in layers:
        layer_bbox = np.asarray(layer.shape.bbox, dtype=float).copy()
        layer_bbox[0] -= bbox_eps
        layer_bbox[1] += bbox_eps
        layer_size = np.maximum(layer_bbox[1] - layer_bbox[0], bbox_eps)
        layer_resolution = np.maximum(
            np.asarray(layer.resolution, dtype=float), bbox_eps
        )
        layer_pattern = patterns.Rectangular.constant(
            element_size=layer_resolution,
            center=layer_bbox[0] + layer_size / 2.0,
        )
        layer_tags = layer_pattern.inside(layer.shape)
        if len(layer_tags) == 0:
            raise RuntimeError(f"No pattern points found inside layer {layer.name!r}.")
        pb.refine_mesh(
            region=layer.shape,
            coordinates=layer_pattern(layer_tags),
        )
        _log(
            f"[build_problem] Layer refinement {layer.name}: "
            f"bbox={layer_bbox.tolist()} res={layer_resolution.tolist()}"
        )

    # Assign materials/boundaries.
    bottom_z = {
        layer.name: float(np.asarray(layer.shape.bbox, dtype=float)[0, 2])
        for layer in layers
    }
    lowest_bottom = min(bottom_z.values())
    region_shapes = {}
    gate_names = []
    for layer in layers:
        _log(
            f"[build_problem] Layer {layer.name}: thickness={layer.thickness}, "
            f"resolution={layer.resolution}, bbox={np.asarray(layer.shape.bbox, dtype=float).tolist()}"
        )

        # Trim shared lower faces off every layer except the lowest one so
        # points that land exactly on an interface have a single owner.
        region = _solver_region(
            layer,
            exclude_bottom=not np.isclose(
                bottom_z[layer.name], lowest_bottom, atol=1e-9, rtol=0.0
            ),
        )
        region_shapes[layer.name] = region

        if _is_gate_layer(layer):
            # Gates are treated as metals with high permittivity.
            pb.set_metal(region=region, setup_name=layer.name)
            pb.set_relative_permittivity(val=1e10, region=region)
            gate_names.append(layer.name)
            _log(f"[build_problem] Gate: {layer.name}")
        elif _is_2deg_layer(layer):
            # 2DEG regions use a selectable boundary condition.
            pb.set_relative_permittivity(val=layer.permitivity, region=region)
            if layer.boundary_condition == "neumann":
                pb.set_neumann(region=region, setup_name=layer.name)
            elif layer.boundary_condition == "helmholtz":
                pb.set_helmholtz(region=region, setup_name=layer.name)
            else:
                pb.set_flexible(region=region, setup_name=layer.name)
            _log(
                f"[build_problem] TwoDEG: {layer.name} "
                f"(bc={layer.boundary_condition}, eps={layer.permitivity})"
            )
        else:
            # Dielectrics use their specified permittivity and a Neumann boundary.
            pb.set_relative_permittivity(val=layer.permitivity, region=region)
            pb.set_neumann(region=region, setup_name=layer.name)
            _log(
                f"[build_problem] Dielectric: {layer.name} "
                f"(eps={layer.permitivity})"
            )

    # Apply Dirichlet boundary on the outer shell of the vacuum.
    inner = shapes.Box(lower_left=corner + vacuum_res, size=size - 2 * vacuum_res)
    boundary = vacuum - inner
    pb.set_dirichlet(region=boundary, setup_name="Boundary")
    region_shapes["Boundary"] = boundary
    _log("[build_problem] Boundary set (Dirichlet).")

    _log("[build_problem] Done.")
    return pb, region_shapes, gate_names

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

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

Parameters:

Name Type Description Default
pescado_problem_finalized Any

Finalized problem instance from Pescado.

 required
gate_shapes RegionMap

Mapping of gate name -> shape.

 required
boundary_shape Shape

Shape for the outer boundary shell.

 required
boundary_name str

Name to use for the boundary in the output mapping.

 'Boundary'

Returns:

Type Description
IndexMap

Dict mapping region names to indices.
Source code in tempura/electrostatics/pescado_wrapper.py 
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
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.

    Args:
        pescado_problem_finalized: Finalized problem instance from Pescado.
        gate_shapes: Mapping of gate name -> shape.
        boundary_shape: Shape for the outer boundary shell.
        boundary_name: Name to use for the boundary in the output mapping.

    Returns:
        Dict mapping region names to indices.
    """
    region_inds: IndexMap = {}
    # "Interior" is everything that is not Dirichlet (Neumann + Helmholtz).
    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

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

Solve the gate basis problem with one MUMPS factorization.

Dirichlet conditions are enforced on the outer boundary (0 V) and gate rows. The raw capacitance matrix is rewritten into a symmetry-preserving Dirichlet system, factored once with MUMPS, and solved in RHS blocks to limit peak memory.

Parameters:

Name Type Description Default
problem_builder ProblemBuilder

Pescado problem builder.

 required
region_shapes RegionMap

Mapping of region names to shapes.

 required
gate_names list[str]

Names of gates to solve.

 required
charge ndarray | None

Optional interior charge columns to include in the RHS.

 None
dtype DTypeLike

Numeric dtype used for factorization and solve.

 DTYPE
blr bool

Whether to enable block low-rank MUMPS factorization.

 BLR
eps_blr float

BLR compression threshold used when blr is enabled.

 EPS_BLR
rhs_block_size int

Number of gate basis columns to solve per MUMPS call.

 8
verbose bool

If True, print progress messages while solving.

 False

Returns:

Type Description
PotentialMap

Tuple of (results, finalized_problem), where results maps gate name to
Any

the solved potential vector.
Source code in tempura/electrostatics/pescado_wrapper.py 
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
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,
    verbose: bool = False,
) -> tuple[PotentialMap, Any]:
    """Solve the gate basis problem with one MUMPS factorization.

    Dirichlet conditions are enforced on the outer boundary (0 V) and gate rows.
    The raw capacitance matrix is rewritten into a symmetry-preserving Dirichlet
    system, factored once with MUMPS, and solved in RHS blocks to limit peak
    memory.

    Args:
        problem_builder: Pescado problem builder.
        region_shapes: Mapping of region names to shapes.
        gate_names: Names of gates to solve.
        charge: Optional interior charge columns to include in the RHS.
        dtype: Numeric dtype used for factorization and solve.
        blr: Whether to enable block low-rank MUMPS factorization.
        eps_blr: BLR compression threshold used when ``blr`` is enabled.
        rhs_block_size: Number of gate basis columns to solve per MUMPS call.
        verbose: If ``True``, print progress messages while solving.

    Returns:
        Tuple of (results, finalized_problem), where results maps gate name to
        the solved potential vector.
    """
    def _log(msg: str) -> None:
        if verbose:
            print(msg)

    _log("[solve_gate_potentials] Finalizing problem...")
    problem = problem_builder.finalized()
    gate_shapes = {name: region_shapes[name] for name in gate_names}
    _log(f"[solve_gate_potentials] Gates: {gate_names}")

    if not gate_names:
        _log("[solve_gate_potentials] No gates provided; returning empty result.")
        return {}, problem

    region_inds = get_region_inds(
        problem,
        gate_shapes,
        region_shapes["Boundary"],
        boundary_name="Boundary",
    )
    boundary_inds = np.asarray(region_inds["Boundary"], dtype=int)
    gate_inds = _ordered_gate_inds(region_inds, gate_names)
    _log(
        f"[solve_gate_potentials] Boundary inds={boundary_inds.size} | "
        f"Gates={len(gate_inds)}"
    )

    dirichlet_inds = _all_dirichlet_inds(boundary_inds, gate_inds)
    A = _symmetric_dirichlet_system(problem.capacitance_matrix, dirichlet_inds)
    _log(
        f"[solve_gate_potentials] System matrix prepared | "
        f"shape={A.shape} | symmetric=True | solver=MUMPS"
    )

    block_size = max(1, min(int(rhs_block_size), len(gate_names)))
    solver = setup_solver(A, dtype=dtype, blr=blr, eps_blr=eps_blr, symmetric=True)
    results: PotentialMap = {}

    # Reuse the factorization across gate blocks; only the prescribed gate
    # columns and optional charge blocks change between solves.
    for start in range(0, len(gate_names), block_size):
        stop = min(start + block_size, len(gate_names))
        gate_block = gate_names[start:stop]
        rhs_input = build_rhs_vectors(
            region_inds,
            gate_block,
            size=A.shape[0],
            charge=_charge_block(
                charge,
                start=start,
                stop=stop,
                total_gates=len(gate_names),
            ),
            dtype=dtype,
        )
        rhs = _dirichlet_rhs(problem.capacitance_matrix, rhs_input, dirichlet_inds)
        _log(
            f"[solve_gate_potentials] Solving block {start}:{stop} | "
            f"rhs_shape={rhs.shape} | rhs_sparse={sp.issparse(rhs)}"
        )

        solutions, solver = solve_linear_system(
            A,
            rhs,
            solver=solver,
            dtype=dtype,
            blr=blr,
            eps_blr=eps_blr,
            symmetric=True,
        )
        if dirichlet_inds.size:
            prescribed = rhs_input[dirichlet_inds, :]
            if sp.issparse(prescribed):
                prescribed = prescribed.toarray()
            solutions[dirichlet_inds, :] = np.asarray(prescribed)

        for i, gate_name in enumerate(gate_block):
            results[gate_name] = np.asarray(solutions[:, i]).copy()

    _log("[solve_gate_potentials] Done.")
    return results, problem
Example