A P-Adaptive Hybridizable Discontinuous Galerkin Spectral Element Method for Electrostatic Particle-in-Cell Simulations

This paper introduces a p-adaptive high-order hybridizable discontinuous Galerkin spectral element method implemented in the open-source PICLas framework to efficiently solve the Poisson equation in electrostatic particle-in-cell simulations by locally refining polynomial degrees in high-gradient regions, thereby significantly reducing the global number of degrees of freedom while effectively modeling complex plasma phenomena.

The Gist

Imagine you are trying to paint a massive, complex mural of a stormy sky. Some parts of the sky are just a smooth, uniform blue (calm weather), while other parts have swirling, chaotic clouds, lightning bolts, and rain (the storm).

The Problem: The "One-Size-Fits-All" Approach
Traditionally, scientists simulating plasma (a super-hot, electrically charged gas found in space thrusters and stars) used a method where they painted the entire mural with the same level of detail. If they wanted to capture the tiny, chaotic lightning bolts accurately, they had to paint the entire sky with that same high level of detail.

This is like using a fine-tipped brush for the whole wall just to get the lightning right. It's incredibly accurate, but it takes forever and uses up all your paint (computing power). Most of the wall (the calm blue sky) doesn't need that much detail, so you are wasting a huge amount of effort.

The Solution: The "Smart Painter" (P-Adaptive HDG-SEM)
This paper introduces a new, smarter way to paint the mural, called a P-Adaptive Hybridizable Discontinuous Galerkin Spectral Element Method. Let's break down what that fancy name means using our painting analogy:

1. The Canvas (The Grid): The computer divides the simulation area into small tiles (like a mosaic).
2. The Brush Strokes (Polynomial Degree): In this method, the "complexity" of the math used to describe the physics on each tile can change.
   • Low Complexity (Low "P"): For the calm blue sky (where the electric field is smooth), the computer uses a simple, low-detail brush. It's fast and cheap.
   • High Complexity (High "P"): For the stormy lightning bolts (where the electric field changes rapidly), the computer automatically switches to a super-fine, high-detail brush.
3. The "Adaptive" Part: The computer acts like a smart painter who constantly looks at the mural. If it sees a storm forming, it instantly upgrades the brush in that specific spot. If the storm clears, it downgrades the brush back to simple. It doesn't waste energy painting the calm sky with a fine brush.

Why is this a big deal?
In the world of space propulsion (like ion thrusters for satellites), we need to simulate how charged particles move. These particles create "sheaths" (layers of charge) that are very thin but have huge changes in electric force.

• Old Way: To see the sheath clearly, you had to make the entire simulation grid very fine. This required millions of calculations, taking days or weeks.
• New Way: The computer only makes the grid fine right where the sheath is. The rest of the simulation stays simple.
  • Result: You get the same accuracy, but you use half the memory and finish the job much faster.

The "Noise" Problem
There was one tricky part. In these simulations, the "paint" is made of individual particles (like millions of tiny dots). Sometimes, just by chance, a few extra dots might cluster together, looking like a storm when it's actually just random noise.
If the computer wasn't smart, it might think, "Oh, look at that cluster! I need a super-fine brush!" and waste time refining a fake storm.
The authors solved this by teaching the computer to distinguish between real physics and statistical noise. They created a rule: "Only upgrade the brush if the pattern is strong enough to be real, not just a random glitch."

The Results
The team tested this "Smart Painter" in three scenarios:

1. A Dielectric Sphere: Like painting a glass ball in a wind. The method correctly used simple brushes inside the glass and complex ones at the surface.
2. A Plasma Sheath: A thin layer of charge near a wall. The method nailed the details of the thin layer while keeping the rest of the simulation fast.
3. An Ion Thruster: A complex 3D engine part. The method successfully zoomed in on the intense electric fields between the grids while ignoring the calmer areas.

In Summary
This paper presents a "smart" simulation tool that knows when to work hard and when to coast. Instead of treating every part of a plasma simulation with the same heavy computational weight, it dynamically adjusts its effort, focusing only where the action is. This allows scientists to simulate complex space engines faster and cheaper, bringing us closer to better space travel technology.

Technical Summary

1. Problem Statement

The simulation of rarefied plasmas, essential for applications like electric space propulsion and semiconductor manufacturing, relies on the Particle-in-Cell (PIC) method to solve the Vlasov-Poisson system. A central computational bottleneck in these simulations is solving the Poisson equation for the electrostatic potential at every time step.

• The Challenge: Traditional PIC codes use low-order finite difference methods. However, plasma phenomena such as sheaths and localized charge distributions create sharp spatial gradients. Resolving these with low-order methods requires extremely fine meshes, leading to prohibitive computational costs and memory usage.
• The Limitation of High-Order Methods: While high-order methods (like Discontinuous Galerkin) offer better accuracy per degree of freedom, applying a uniform high polynomial degree ($p$) across the entire domain is inefficient because high resolution is only needed in specific regions with strong gradients.

2. Methodology

The authors propose a $p$-adaptive Hybridizable Discontinuous Galerkin Spectral Element Method (HDG-SEM) integrated into the open-source framework PICLas.

A. Numerical Formulation (HDG-SEM)

• Hybridization: The method introduces a trace variable ($\lambda$) representing the electric potential on element interfaces. This allows for the static condensation of interior degrees of freedom, resulting in a global system that only solves for the interface potentials.
• Spectral Elements: The domain is discretized into hexahedral elements using tensor-product polynomial bases (Lagrange polynomials with Legendre-Gauss nodes).
• Weak Formulation: The Poisson equation is transformed into a first-order system involving the electric displacement field ($D$) and potential ($\phi$). Numerical fluxes are defined using a stabilization parameter ($\tau$) to handle discontinuities at element interfaces.
• Integration: Gaussian quadrature is used for integration. The method employs a bijective mapping from reference elements to physical elements to handle curved boundaries.

B. P-Adaptation Strategy

The core innovation is the automatic, element-local adjustment of the polynomial degree ($k$) based on the solution's local smoothness.

• Modal Analysis: The solution (electric potential $\phi$ and charge densities $\rho$) within each element is expanded into an orthogonal Legendre basis.
• Refinement Indicator:
  • For smooth fields (potential), the polynomial degree is increased if the modal coefficients decay slowly.
  • For charge densities (derived from discrete particles), statistical noise is a major issue. The authors introduce a dynamic threshold based on the estimated variance ($\sigma^2_k$) of the modal coefficients caused by particle sampling.
  • Refinement only occurs if the physical signal significantly exceeds the statistical noise ($\hat{\rho}_k^2 > \epsilon \sigma_k^2$). This prevents the solver from wasting resources resolving particle noise.
• Implementation: The data structure uses derived types to handle variable polynomial degrees per element and face. The global system is solved using PETSc (with Cholesky decomposition for small systems or Conjugate Gradient with BoomerAMG preconditioning for large systems).

3. Key Contributions

1. Novel Algorithm: Development of a $p$-adaptive HDG-SEM specifically tailored for electrostatic PIC simulations, addressing the unique challenge of statistical noise in particle-based charge deposition.
2. Noise-Aware Adaptation: Introduction of a variance-based thresholding mechanism that distinguishes between physical gradients and statistical particle noise, preventing spurious refinement.
3. Efficiency Gains: Demonstration that the method can concentrate computational effort strictly on regions with high gradients (e.g., plasma sheaths) while using low-order approximations in smooth regions, significantly reducing the global degrees of freedom (DOFs).
4. Open-Source Integration: Successful implementation and validation within the massively parallel PICLas framework.

4. Validation Results

The method was validated against three benchmark cases:

• Dielectric Sphere (Analytical Benchmark):

  • Setup: A sphere in a uniform electric field.
  • Result: The $p$-adaptive solver achieved accuracy comparable to uniform high-order simulations but used a lower polynomial degree ($N=2$) inside the sphere where the field is linear, and higher degrees only near the interface. This reduced the total DOFs significantly.

• 1D Plasma Sheath (Coupled PIC):

  • Setup: A collisionless plasma interacting with a charged wall, creating a sharp potential gradient.
  • Comparison:
    • Uniform Low-Order (N=1): Failed to resolve the sheath.
    • Uniform High-Order (N=4): Accurate but used 20 DOFs.
    • P-Adaptive: Used $N=4$ only in the sheath region and $N=1$ elsewhere.
  • Result: The adaptive method achieved similar accuracy to the uniform high-order case with 14 DOFs (vs. 20), demonstrating a ~30% reduction in computational cost for the same accuracy. It outperformed $h$-refinement (increasing element count), which required 64 DOFs for similar accuracy.

• 2D Axisymmetric Ion Optic (Complex Geometry):

  • Setup: Simulation of a gridded ion thruster aperture with strong localized fields between grids.
  • Result: The solver successfully identified the region between the grids (where gradients are high) and increased the polynomial degree up to $N=5$. In the bulk plasma and downstream regions, it maintained $N=1$.
  • Limitation Observed: In regions with very low particle counts (near the symmetry axis), the variance estimate was high, preventing refinement despite physical gradients. This highlights a trade-off between noise filtering and resolution in low-density regions.

5. Significance and Future Work

• Significance: This work bridges the gap between high-order accuracy and computational efficiency in PIC simulations. By dynamically adapting the polynomial order, it enables the simulation of complex plasma phenomena with fewer resources, making high-fidelity simulations of space propulsion and industrial plasma processes more feasible.
• Future Directions:
  • Developing a fully dynamic adaptation strategy that adjusts $p$ during transient simulations (currently done iteratively between runs).
  • Refining the variance estimation to account for non-uniform particle distributions and residence times.
  • Extending the method to 3D unstructured meshes with complex chemical reactions and ionization.

In conclusion, the paper presents a robust, noise-aware, $p$-adaptive solver that significantly optimizes the computational cost of electrostatic PIC simulations without sacrificing accuracy in critical gradient regions.