Surface differential operators on the torus¶

Two models in this repository require surface differential operators:

The electrode model solves a Poisson equation on the surface (Laplace–Beltrami operator).
The current-potential model computes \(\nabla_s \Phi\) to obtain a divergence-free surface current.

Surface gradient¶

For a scalar field \(f(\theta,\phi)\) on a parameterized surface, the surface gradient can be written as:

\[ \nabla_s f = g^{ij} \frac{\partial f}{\partial u^i}\, \mathbf r_{u^j}, \]

where \(u^1=\theta\), \(u^2=\phi\), and \(g^{ij}\) is the inverse metric.

For the circular torus (orthogonal coordinates, \(F=0\)):

\[ \nabla_s f = \frac{1}{E}\, f_\theta\, \mathbf r_\theta + \frac{1}{G}\, f_\phi\, \mathbf r_\phi = \frac{1}{a^2}\, f_\theta\, \mathbf r_\theta + \frac{1}{R(\theta)^2}\, f_\phi\, \mathbf r_\phi. \]

Laplace–Beltrami (surface Laplacian)¶

The Laplace–Beltrami operator is:

\[ \Delta_s f = \nabla_s \cdot \nabla_s f = \frac{1}{\sqrt{g}} \left[ \frac{\partial}{\partial \theta} \left(\sqrt{g}\, g^{\theta\theta}\, f_\theta\right) + \frac{\partial}{\partial \phi} \left(\sqrt{g}\, g^{\phi\phi}\, f_\phi\right) \right]. \]

For the circular torus, with \(\sqrt{g}=aR(\theta)\), \(g^{\theta\theta}=1/a^2\), \(g^{\phi\phi}=1/R(\theta)^2\), this becomes:

\[ \Delta_s f = \frac{1}{a^2 R(\theta)} \frac{\partial}{\partial \theta}\left(R(\theta)\, f_\theta\right) + \frac{1}{R(\theta)^2}\, f_{\phi\phi}. \]

Discretization and numerics in this code¶

Periodic grids¶

We discretize \(\theta\) and \(\phi\) on uniform periodic grids:

\[ \theta_j = \frac{2\pi j}{N_\theta},\qquad \phi_k = \frac{2\pi k}{N_\phi}. \]

Spectral derivatives¶

Derivatives in \(\theta\) and \(\phi\) are computed using FFTs:

\[ \frac{\partial f}{\partial x} \approx \mathcal{F}^{-1}\left( i k\, \mathcal{F}(f)\right), \]

implemented in src/torus_solver/spectral.py.

Poisson solve¶

Because \(R(\theta)\) varies, \(\Delta_s\) is not diagonal in Fourier space, so the Poisson equation is solved with Conjugate Gradient (CG) using a matrix-free operator.

In the electrode model we solve (schematically):

\[ -\Delta_s V = s/\sigma_s, \]

with:

a net-zero constraint on \(s\) (closed surface)
a gauge fix on \(V\) (remove the constant nullspace)

See src/torus_solver/poisson.py.