Electrostatics in 2D

Setup

The simulation solves for the electric potential \(V\) and electric field \(\mathbf{E} = -\mathbf{\nabla} V\) in a 2D region with various charge distributions and dielectric materials and ideal conductors. The potential is determined by Gauss' law $$ \mathbf{\nabla} \cdot \mathbf{D} = \rho_f $$ where \(\mathbf{D} = \epsilon \mathbf{E}\) is the electric displacement field, \(\epsilon\) is the permittivity of the material, and \(\rho\) is the free charge density. Ideal conductors are treated as regions of fixed potential. We work in units where the permittivity of free space is unity, \(\epsilon_0 = 1\). In terms of the potential this reads $$ -\mathbf{\nabla} \cdot (\epsilon \mathbf{\nabla} V) = \rho_f $$ In two dimensions this can be expanded as $$ \frac{\partial}{\partial x} \left( \epsilon \frac{\partial V}{\partial x} \right) + \frac{\partial}{\partial y} \left( \epsilon \frac{\partial V}{\partial y} \right) = -\rho_f $$ The electric field is related to the potential by \(\mathbf{E} = -\mathbf{\nabla} V\) and the charge density (free and bound) is determined from the electric field as \(\rho = \rho_f + \rho_b = \mathbf{\nabla} \cdot \mathbf{E}\). The outer boundaries of the simulation are grounded (\(V=0\)).

Method

We discretize these equations on a square grid \(V(\mathbf{r}) \rightarrow V_{i,j}\) and use finite differences to approximate the derivatives. This yields a set of linear equations for the \(V_{i,j}\). These linear equations are solved iteratively using a Gauss-Seidel method with successive over-relaxation (SOR). To parallelize the computation on the GPU, we use a two pass red-black scheme. The over-relaxation parameter is set to \(\omega = 1.985\).