Multigrid Poisson Solver for Complex Geometries Using Finite Difference Method

We present an efficient numerical method, inspired by transformation optics, for solving the Poisson equation in complex and arbitrarily shaped geometries. The approach operates by mapping the physical domain to a uniform computational domain through coordinate transformations, which can be applied either to the entire domain or selectively to specific boundaries inside the domain. This flexibility allows both homogeneous (Laplace equation) and inhomogeneous (Poisson equation) problems to be solved efficiently using iterative or fast direct solvers, with only the material parameters and source terms modified according to the transformation. The method is formulated within a finite difference framework, where the modified material properties are computed from the coordinate transformation equations, either analytically or numerically. This enables accurate treatment of arbitrary geometric shapes while retaining the simplicity of a uniform grid solver. Numerical experiments confirm that the method achieves second-order accuracy , and offers a straightforward pathway to integrate fast solvers such as multigrid methods on the uniform computational grid.

💡 Research Summary

The paper introduces a novel numerical framework for solving the Poisson equation in domains with complex, arbitrarily shaped boundaries. Inspired by transformation optics, the authors map the physical (non‑uniform) domain onto a uniform computational domain through a coordinate transformation x → x′. The Jacobian J  of this transformation and its determinant are used to construct a spatially varying permittivity tensor ε in the computational space, while the physical domain is assigned a unit tensor (ε′ = I). This mapping preserves the original physics: the Poisson equation ∇·(ε∇φ)=−ρ in the computational domain is mathematically equivalent to the original equation in the physical domain, with the source term scaled by det(J).

The formulation is presented in two dimensions for clarity, but the authors note that extension to three dimensions is straightforward. Explicit expressions for the components of ε (ε_xx, ε_xy, ε_yx, ε_yy) are derived from the Jacobian entries. The electric field is expressed as the gradient of the potential (E = −∇φ), leading to a generalized Poisson equation that incorporates the anisotropic ε. A second‑order central finite‑difference stencil is then applied, yielding an update formula for the potential at each grid point that involves six coefficients (k₁–k₆) dependent on ε and the grid spacings Δx, Δy. The relaxation factor ω (1 < ω < 2) provides a Gauss‑Seidel smoothing step; ω = 1 reduces to the classic Gauss‑Seidel method, which is used as the smoother within the multigrid cycle.

The core of the contribution lies in integrating this transformed‑domain finite‑difference discretization with a geometric multigrid (GMG) solver. Because the computational grid is uniform, the standard multigrid hierarchy (fine, coarse, coarser grids) can be generated by simple halving of the mesh spacing, and restriction/prolongation operators are the conventional full‑weighting and bilinear interpolation. The residual equation A_h Ψ_h = r_h is solved on progressively coarser levels, and the correction is interpolated back to the fine grid. Since ε and the source term have already been incorporated into the system matrix A_h, no additional matrix assembly or complex coarsening of an unstructured mesh is required. This results in an O(N) algorithm with minimal preprocessing overhead and excellent data locality, making it well‑suited for parallel implementations on modern CPUs or GPUs.

To validate the method, two sets of numerical experiments are presented. The first uses an analytically defined mapping that transforms a rectangular computational domain into a semi‑annular physical domain. The mapping x′ = x cos y, y′ = x sin y produces a non‑uniform ε with only diagonal components (ε_xy = ε_yx = 0). The analytical solution for the potential in this geometry is known (logarithmic radial dependence). Simulations on a 201 × 201 grid demonstrate second‑order convergence in L∞, L₁, and L₂ norms, with the error curves aligning with the Δx² trend. Moreover, the multigrid solver achieves a speed‑up of roughly two orders of magnitude compared to a Successive Over‑Relaxation (SOR) solver (ω = 1.875, tolerance = 10⁻¹⁰) on the same hardware (Intel Core i7‑1165G7, 16 GB RAM).

The second set addresses cases where no closed‑form mapping is available. Here, the physical grid is generated numerically (e.g., via elliptic grid generation or other mesh‑generation tools), and the Jacobian is computed numerically using finite differences. The same workflow—construction of ε, finite‑difference discretization, multigrid solution, and back‑transformation of the potential—remains applicable. The authors report that the method retains its second‑order accuracy and comparable performance gains, confirming its robustness for arbitrary geometries.

In the discussion, the authors highlight several advantages: (1) the ability to reuse existing structured‑grid FD and multigrid codes without substantial modification; (2) preservation of physical quantities (potential, electric field) through exact transformation formulas; (3) optimal computational complexity (O(N)) and high parallel efficiency; and (4) avoidance of the assembly and storage costs associated with FEM on unstructured meshes. Potential limitations are also acknowledged: highly distorted Jacobians can lead to large variations in ε, which may cause numerical stiffness or require additional smoothing/preconditioning strategies. Nonetheless, the overall conclusion is that transformation‑optics‑based coordinate mapping provides a practical, accurate, and fast alternative to traditional FEM‑based solvers for Poisson problems in complex domains.