Scribd
Upload
5 views3 pages

2D Electrostatic PIC Algorithm Guide

The document outlines a generic 2D electrostatic/electromagnetic Particle-In-Cell (PIC) algorithm, detailing initialization, main time-step loop, and post-processing steps. Key components include particle injection, charge/current deposition, field solving, particle pushing, and diagnostics. It emphasizes the importance of self-consistent fields, resolution conditions, and energy conservation in the simulation process.

Uploaded by

Md Arifuzzaman Faisal
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
5 views3 pages

2D Electrostatic PIC Algorithm Guide

The document outlines a generic 2D electrostatic/electromagnetic Particle-In-Cell (PIC) algorithm, detailing initialization, main time-step loop, and post-processing steps. Key components include particle injection, charge/current deposition, field solving, particle pushing, and diagnostics. It emphasizes the importance of self-consistent fields, resolution conditions, and energy conservation in the simulation process.

Uploaded by

Md Arifuzzaman Faisal
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd

PIC Algorithm

Here’s a generic 2D electrostatic / electromagnetic PIC algorithm, laid out in logical steps and
pseudocode. You can adapt it to your FDTD–PIC coupling for Smith–Purcell:

1. Initialization
1. Define physical constants
1
ε 0 , μ0 , c= , me , q e , …
√ ε0 μ 0
2. Set up spatial grid

o $N_{x},\text{\:\,}N_{y}$ cells, spacing $\Delta x,\text{\:\,}\Delta y$


o Field arrays: E ( i , j ), B (i , j ) (staggered Yee grid)
o Charge density array ρ ( i , j ), current density J ( i, j )
3. Load or define conductor/geometry masks (e.g. grating)

4. Initialize particles

o Positions ( x p , y p ), velocities ( v x , p , v y, p , v z , p )
o Flags “active/inactive”
5. Precompute any static fields (e.g. Bx ) and PML damping profiles

2. Main time-step loop (for n=1 … N s t e p s)

2.1 Particle injection (if continuous beam)


for each injection slot:
find inactive particle p
x_p ← injection plane + small offset
y_p ← beam center + random spread
v_{x,p} ← beam speed, v_{y,p}=0, v_{z,p}=0
mark p active
end

2.2 Charge/current deposition (“scatter”)


1. Zero ρ ( i , j ), J ( i, j )

2. For each active particle p:

o Locate cell indices ( i , j ) ≃ $\left\lfloor (x_{p} - x_{\min})/\Delta x \


right\rfloor,\text{\:\,}\lfloor(y_{p} - y_{\min})/\Delta y\rfloor$
o Compute weights W i ± 1, j ± 1 (e.g. linear or quadratic)

o Accumulate

$$\rho(i,j)\text{\:\,} + \text{ ⁣} = \text{\:\,}q_{e}\, W_{i,j},\quad J_{x}(i,j)\text{\:\,} + \text{ ⁣} = \


text{\:\,}q_{e}\, v_{x,p}\, W_{i,j},\text{\:\,}\ldots$$

2.3 Field solve (FDTD or Poisson)


 Electromagnetic (Yee FDTD):
B^{n+½} ← B^{n−½} − Δt ∇×E^n
E^{n+1} ← E^n + Δt ( c² ∇×B^{n+½} − J/ε₀ )

 Electrostatic (if applicable): solve ∇²φ=−ρ/ε₀ then E=−∇φ

 Apply PML or boundary conditions

2.4 Field gather to particles (“interpolate”)


For each active p:
1. Locate surrounding grid nodes ( i , j )
2. Interpolate E p =∑ E ( i, j ) W i , j, B p=∑ B ( i, j ) W i , j

2.5 Particle push (Boris scheme)


v⁻ = v^n + (q/m) E_p (Δt/2)
t = (q/m) B_p (Δt/2)
s = 2 t / (1+|t|²)
v′ = v⁻ + v⁻×t
v⁺ = v⁻ + v′×s
v^{n+1} = v⁺ + (q/m) E_p (Δt/2)
x^{n+1}_p = x^n_p + v^{n+1} Δt

 Handle domain exits → deactivate


2.6 Diagnostics & visualization
 Track energy conservation, particle counts
 Plot fields/particles every few steps

3. End of loop / post-processing


 Compute spectra of fields (e.g., via FFT)
 Extract radiation intensity vs. angle or frequency

Key points
 Charge → field coupling via deposition and interpolation ensures self-consistent fields.
 Lambda_D (Debye length) and c/ωₚ must be resolved by Δ x , Δ t .
 Courant condition: Δ t <1/ ( c √ 1/ Δ x 2+ 1/ Δ y 2 ).
 Boris push is symplectic and preserves energy very well.
 PML absorbs outgoing waves to avoid reflections.
Feel free to adapt each step to your vectorized MATLAB implementation and to your metallic‐
grating boundary treatment.

Common questions

Powered by AI

The process of field gathering involves locating the grid nodes surrounding each active particle and interpolating the electric and magnetic field values from these nodes onto the particle's position. This interpolation uses weights calculated for the surrounding grid positions, ensuring that each particle interaction with the field is consistent with its location relative to the grid .

In the 2D PIC algorithm, charge and current deposition is handled by first zeroing out the charge density ρ(i, j) and current density J(i, j) arrays. Then, for each active particle, the indices of the cell it resides within are located, and weights are computed for neighboring cells. These weights are used to accumulate charge and current densities in the respective cell locations .

The Boris scheme is symplectic, meaning it preserves the phase space volume and energy very well by effectively splitting the updates into half time-step electric field integrations, rotating via a magnetic field, and completing the full time-step update with another half time-step electric field integration, thus ensuring energy is conserved over long simulation times .

Perfectly Matched Layers are implemented to absorb outgoing waves at the boundaries of the simulation domain. By absorbing these waves, PMLs prevent reflection back into the computational domain, which can cause errors and artificial interference in the simulation, improving the accuracy of the results in open-boundary scenarios .

Setting up conductor/geometry masks is crucial in defining the spatial distribution and effects of boundaries within the simulation. For a metallic grating, these masks determine how electromagnetic fields interact with physical obstacles, impacting scattering, absorption, and other boundary-related phenomena. This allows the simulation to accurately reflect real-world scenarios involving complex geometries .

Diagnostics in PIC simulations are used to track energy conservation and particle count, provide visualizations of fields and particles at intervals, and help in assessing the accuracy of simulation by comparing expected results. Monitoring these diagnostics helps in adjusting simulation parameters and methods to ensure that the simulation remains accurate and stable over time .

The Courant condition imposes a stability requirement on the timestep, stating that Δt must be less than 1/(c√1/Δx² + 1/Δy²), where c is the speed of light. This ensures that numerical wave propagation does not exceed the physical wave speed, preventing instabilities and ensuring accurate time evolution of the fields .

Resolving the Debye length (λ₄) and the plasma frequency (ωₚ) ensures that the grid spacing (Δx) and timestep (Δt) are sufficiently small to capture the physical phenomena correctly. The Debye length is crucial for resolving the shortest scale of charge separation, and the plasma frequency determines the fastest natural oscillations. Failing to resolve these scales can lead to numerical instability and inaccurate representation of physical behaviors .

The necessary initializations include defining physical constants such as ε₀, μ₀, and c, setting up the spatial grid with Nx, Ny cells and spacing Δx, Δy, initializing field arrays E(i, j), B(i, j) on a staggered Yee grid, and setting up charge density array ρ(i, j) and current density J(i, j). Additionally, particles should be initialized with positions, velocities, and active/inactive flags, and any static fields and PML damping profiles should be precomputed .

Within the PIC algorithm framework, electromagnetic field solutions often use the Yee FDTD method to update magnetic and electric fields in time using curl operations and current density contributions. Electrostatic solutions, on the other hand, solve Poisson's equation for the potential, which is then used to derive electric fields by calculating the negative gradient of the potential. The choice between these methods typically depends on relative speed and field configurations in the simulation .

You might also like