This page summarizes the mathematical model used by Tempura and the solver assumptions behind the current Poisson workflow.

Continuum model

The electrostatic potential \(\phi(\mathbf{r})\) is modeled by the Poisson equation

\[ -\nabla \cdot \bigl(\varepsilon(\mathbf{r}) \, \nabla \phi(\mathbf{r})\bigr) = \rho(\mathbf{r}), \]

where:

  • \(\phi\) is the electrostatic potential
  • \(\varepsilon\) is the spatially varying permittivity
  • \(\rho\) is the charge density

In Tempura, the layer stack defines the geometry and the material map \(\varepsilon(\mathbf{r})\). Pescado then meshes that geometry and assembles the discrete sparse linear system.

From geometry to a sparse linear system

After the device is finalized and meshed, the electrostatic problem is reduced to a linear system

\[ A x = b, \]

where:

  • \(A\) is the sparse matrix assembled by Pescado
  • \(x\) is the vector of unknown electrostatic potentials on the mesh
  • \(b\) collects source terms and prescribed-voltage contributions

In the current implementation, the assembled matrix is exposed as problem.capacitance_matrix. The gate-basis workflow keeps the geometry, material map, and response law fixed, and changes only the prescribed gate voltages.

Boundary-condition roles in Tempura

Tempura uses four boundary-condition classes during problem assembly.

Dirichlet

Dirichlet conditions prescribe the potential directly:

\[ \phi = \phi_0. \]

In this codebase, Dirichlet constraints are used for:

  • the outer vacuum shell, which is grounded
  • gate regions, which are prescribed to fixed voltages during basis solves

Neumann

Neumann conditions prescribe the normal flux:

\[ \hat{\mathbf{n}} \cdot \bigl(\varepsilon \, \nabla \phi\bigr) = g. \]

In Tempura, dielectric regions are typically tagged with Neumann conditions.

Helmholtz

The Helmholtz condition used by Pescado represents a linear response law, often written as

\[ \hat{\mathbf{n}} \cdot \bigl(\varepsilon \, \nabla \phi\bigr) + \alpha \, \phi = \beta. \]

This is the default model for the TwoDEG layer because it preserves a linear relationship between potential and response.

Flexible

Flexible regions are placeholders that can later be reassigned to Dirichlet, Neumann, or Helmholtz subsets after meshing.

Gate-basis solve

solve_gate_potentials(...) computes one unit-voltage basis solution per gate. For a chosen gate \(i\), the solve sets:

  • gate \(i\) to \(1\,\mathrm{V}\)
  • all other gates to \(0\,\mathrm{V}\)
  • the outer boundary to \(0\,\mathrm{V}\)

If the model remains linear, any gate-voltage configuration can then be formed by linear superposition of those basis solutions.

This is why the current workflow is efficient: the expensive geometry assembly and matrix factorization are done once, and each additional gate only adds another right-hand side.

Symmetric Dirichlet rewrite

The gate and boundary nodes are prescribed values, so they are not free unknowns in the usual sense. Let the unknowns be split into:

  • interior indices \(I\)
  • constrained Dirichlet indices \(D\)

Then the original system has the block form

\[ \begin{bmatrix} A_{II} & A_{ID} \\ A_{DI} & A_{DD} \end{bmatrix} \begin{bmatrix} x_I \\ x_D \end{bmatrix} = \begin{bmatrix} b_I \\ b_D \end{bmatrix}, \qquad x_D = g. \]

The interior equation is therefore

\[ A_{II} x_I = b_I - A_{ID} g. \]

Tempura does not physically remove the constrained rows and columns from the matrix. Instead, it rewrites the full system into an equivalent constrained system that preserves symmetry:

  • Dirichlet rows are zeroed
  • Dirichlet columns are zeroed
  • Dirichlet diagonal entries are set to 1
  • the removed coupling term \(A_{:D} g\) is shifted to the right-hand side

This is the role of _symmetric_dirichlet_system(...) and _dirichlet_rhs(...) in the solver module.

Why preserve symmetry? Because the current MUMPS path is called with symmetric=True. If only the Dirichlet rows were replaced and the columns were left untouched, the modified matrix would generally become non-symmetric and would no longer match the factorization mode.

Role of MUMPS

MUMPS is the sparse direct linear solver used for the finalized constrained system.

In the current implementation, MUMPS is used in three steps:

  1. build the constrained sparse matrix once
  2. factor that matrix once
  3. solve for many right-hand sides using the same factorization

This is a good fit for the gate-basis workflow because all basis solves share the same matrix and differ only in the prescribed gate columns and optional charge columns.

Tempura also exposes optional block low-rank (BLR) settings when MUMPS is used. Those settings affect factorization performance and memory use, but they do not change the mathematical problem being solved.

Matrix constraints in the current solver path

The current implementation assumes the following:

  • the assembled matrix is square
  • the matrix is sparse enough to benefit from sparse direct factorization
  • after Dirichlet rewriting, the matrix remains symmetric because the solver is invoked with symmetric=True
  • the set of Dirichlet indices is consistent across the solve and corresponds to the outer boundary plus the gate nodes
  • any supplied charge array has the same number of rows as the linear system, and either one column or one column per gate block

There are also model-level constraints:

  • linear superposition is valid only while the geometry, material map, and linearized response law remain fixed
  • if a future model introduces nonlinear charge response, voltage-dependent material properties, or a non-symmetric discretization, then the current "factor once, solve many RHS" approach will no longer apply directly

Practical interpretation

In practical terms, Tempura solves Poisson's equation by:

  1. converting the layer stack into a meshed electrostatic domain with Pescado
  2. assembling one sparse matrix for that fixed geometry
  3. rewriting the matrix to enforce Dirichlet gate and boundary constraints while preserving symmetry
  4. factoring the constrained matrix once with MUMPS
  5. solving one right-hand side per gate basis vector, possibly in blocks

That is the implementation-level meaning of the current Poisson solve in Tempura.