Particle Tracking and Charge Accumulation

7. Particle generation and injection state

Source: init_particle_batch_from_config, bem_injection, bem_sheath_runtime

7.1 Species processing

init_particle_batch_from_config traverses all species, determines the number generated on each rank, and interleaves them into the SoA particle array. In MPI runs, mpi_split_count distributes the global count across ranks.

Supported source_mode values:

source_mode	Count	Position	Velocity
volume_seed	npcls_per_step	Uniform random in pos_low to pos_high	Shifted Maxwellian
reservoir_face	Determined dynamically from flux and batch duration	Injection-face rectangle	Maxwellian weighted by incoming flux, or velocity grid
photo_raycast	Number of hits from rays_per_batch rays	First surface hit by each ray	Flux-weighted Maxwellian along the surface normal

7.2 reservoir_face macro count

The incoming flux of a drifting Maxwellian is computed from the normal velocity component against the injection-face inward normal n:

$$u_n = \mathbf{u} \cdot \mathbf{n}$$

and the thermal speed:

$$\sigma = \sqrt{\frac{k_B T}{m}}$$

The implementation uses flux_weighted_normal_tail(vmin, u_n, sigma) and counts only particles whose normal velocity is at least vmin_normal.

For area A, batch duration T_b, and macro-particle weight w:

$$N_\mathrm{phys} = \Gamma_\mathrm{in} A T_b$$ $$N_\mathrm{macro,expected} = \frac{N_\mathrm{phys}}{w}$$

The fractional part is carried in injection_state%macro_residual(species).

$$B = r_\mathrm{old} + N_\mathrm{macro,expected}$$ $$N_\mathrm{macro} = \lfloor B \rfloor,\quad r_\mathrm{new} = B - \lfloor B \rfloor$$

In MPI runs, batch_duration_scale = 1 / nrank, so each rank generates its share of the global flux.

7.3 Reservoir sampling

For reservoir_face, positions are sampled uniformly from the injection-face rectangle and, if requested, shifted slightly away from the face by position_jitter_dt=sim.dt. Velocity is generated by one of the following methods:

• Generate a shifted Maxwellian and resample only the normal component from a flux-weighted distribution.
• If velocity_distribution="grid", read a CSV velocity grid. With phase_space, use max(v_n,0) f(v); with flux_weighted, treat the input values as the incoming distribution.

With reservoir_potential_model="infinity_barrier", the normal-velocity lower bound is corrected from the difference between the injection-face mean potential and phi_infty.

7.4 photo_raycast

photo_raycast launches rays from the injection face and emits particles from the first mesh element hit by each ray.

1. Sample the ray origin uniformly on the injection-face rectangle.
2. Normalize ray_direction and confirm that it points inward.
3. Extend the ray to the box boundary and run find_first_hit on that segment.
4. If it hits the mesh, use the element normal to construct emission position and velocity.
5. If it hits a box boundary and the boundary condition reflects or periodically wraps it, continue up to raycast_max_bounce.

The weight per hit is:

$$w_\mathrm{hit} = \frac{J_\perp A_\perp T_b} {|q| N_\mathrm{ray,global}}$$

where A_perp = A * abs(dot(ray_direction, inward_normal)).

If deposit_opposite_charge_on_emit=true, the source element receives:

$$\Delta q_\mathrm{emit} = -q_\mathrm{particle} w_\mathrm{eff}$$

as photo_emission_dq.

7.5 Photo escape closure

With photo_escape_model="boltzmann_cutoff", BEACH uses the centroid potential excluding the contribution from the emitting element itself.

$$\mathrm{barrier} = \max(\phi_\mathrm{emit} - \phi_\infty, 0)$$ $$f_\mathrm{escape} = \exp\left( {}-\frac{|q|\,\mathrm{barrier}}{k_B T_\mathrm{PE}} \right)$$

The effective weight is:

$$w_\mathrm{eff} = w_\mathrm{hit} f_\mathrm{escape}$$

This is a reduced closure: returning photoelectrons are not tracked individually and are treated as immediate neutralization.

---

8. Boris pusher

Source: bem_pusher

Particle motion is advanced with the Boris method using the uniform magnetic field sim.b0 and the electric field E evaluated at the particle position.

Inputs:

• position x
• velocity v
• charge q
• mass m
• time step dt
• electric field E
• magnetic flux density B

Update equations:

$$\mathbf{v}^- = \mathbf{v}^n {}+ \frac{q}{m}\mathbf{E}\frac{\Delta t}{2}$$ $$\mathbf{t} = \frac{q}{m}\mathbf{B}\frac{\Delta t}{2} ,\quad \mathbf{s} = \frac{2\mathbf{t}}{1 + \lVert\mathbf{t}\rVert^2}$$ $$\mathbf{v}' = \mathbf{v}^- + \mathbf{v}^- \times \mathbf{t}$$ $$\mathbf{v}^+ = \mathbf{v}^- + \mathbf{v}' \times \mathbf{s}$$ $$\mathbf{v}^{n+1} = \mathbf{v}^+ + \frac{q}{m}\mathbf{E}\frac{\Delta t}{2}$$ $$\mathbf{x}^{n+1} = \mathbf{x}^{n} + \mathbf{v}^{n+1}\Delta t$$

BEACH runs collision detection on the segment x^n -> x^{n+1}. If there is a collision, the particle is absorbed and x^{n+1} is not saved back into particle state.

---

9. Collision detection

Source: bem_collision, bem_mesh

9.1 Broad phase

init_mesh builds each element AABB and the collision grid. Small meshes use linear search; larger meshes use a uniform grid plus 3D-DDA.

Collision grid:

1. Build the bounding box of all element AABBs.
2. Estimate cell width from target_elems_per_cell.
3. Register each element index, in CSR form, into cells overlapped by its AABB.

The particle segment p0 -> p1 is first intersected with the grid AABB, and only traversed cells are enumerated by 3D-DDA. Only elements registered in those cells are passed to the narrow phase.

9.2 Narrow phase: Moller-Trumbore

The segment

$$\mathbf{p}(t) = \mathbf{p}_0 + t(\mathbf{p}_1 - \mathbf{p}_0), \quad 0 \le t \le 1$$

and the triangle

$$\mathbf{v}(u,v) = \mathbf{v}_0 + u(\mathbf{v}_1-\mathbf{v}_0) + v(\mathbf{v}_2-\mathbf{v}_0)$$

are tested with the Moller-Trumbore method. The conditions are:

• the triangle is not degenerate
• the segment direction is not nearly parallel to the triangle plane
• 0 <= u <= 1
• 0 <= v
• u + v <= 1
• 0 <= t <= 1

If multiple elements are hit, the one with the smallest t is selected.

9.3 periodic2 collision

In periodic2, the mesh itself contains only base elements. find_first_hit_periodic2 computes the required image-shift range from the segment and mesh AABB.

$$n_\mathrm{min} = \left\lceil \frac{\operatorname{min}(p_0, p_1) - \operatorname{max}(\mathrm{mesh}) - \mathrm{tol}}{L} \right\rceil$$ $$n_\mathrm{max} = \left\lfloor \frac{\operatorname{max}(p_0, p_1) - \operatorname{min}(\mathrm{mesh}) + \mathrm{tol}}{L} \right\rfloor$$

For each image, the segment is shifted by -shift and intersected against the base mesh. The hit record stores both the physical image-coordinate hit position hit%pos and the primary-cell wrapped position hit%pos_wrapped.

---

10. Box boundary conditions

Source: bem_boundary

If a particle does not hit the mesh and the candidate updated position leaves the simulation box, axis-wise boundary conditions are applied.

Boundary condition	Processing
open	Remove the particle and count it as escaped_boundary
reflect	Mirror the position at the boundary plane and reverse the normal velocity component
periodic	Wrap to the opposite side

apply_box_boundary checks the three axes in order. If one step crosses multiple periodic lengths, periodic still returns the particle to the box using modulo. For reflect and periodic, the result is clamped slightly inside the box, around 1e-12, to avoid numerical instability from landing exactly on a boundary.

---

11. Charge deposition and surface models

Source: commit_batch_charge, bem_surface_models

11.1 Insulator accumulation

With the default surface_model="insulator", absorbed particle charge is accumulated directly on the hit element.

When particle p hits element i:

$$\Delta q_i \mathrel{+}= q_p w_p$$

This accumulated charge is used as source charge in the next batch’s field-solver refresh.

11.2 Photo emission bookkeeping

When photo_raycast uses deposit_opposite_charge_on_emit=true, the opposite-sign charge is added to the emitting element. This is accumulated as photo_emission_dq, separately from charge deposited by later particle collisions, and merged at batch commit.

11.3 Floating conductor relaxation

surface_model="conductor" is available only when field_bc_mode="free". Conductor elements are grouped by mesh_id as floating conductor groups. The goal is to equalize element potentials within each conductor object while conserving the total charge of that object.

Unknowns:

• conductor element charges q_i
• equipotential value V_g for each conductor group

For an element i in group g(i):

$$\sum_j A_{ij} q_j - V_{g(i)} = {}-\phi_i^\mathrm{fixed}$$

where

$$A_{ij} = \begin{cases} 1/\epsilon, & i=j,\ \epsilon>0, \\ 2\sqrt{\pi}/h_i, & i=j,\ \epsilon=0, \\ 1/\sqrt{\lVert\mathbf{c}_i-\mathbf{c}_j\rVert^2+\epsilon^2}, & i\ne j \end{cases}$$

phi_fixed is the potential from non-conductor charge and the uniform external electric field, divided by k_coulomb.

Each group adds a total-charge conservation constraint:

$$\sum_{i \in g} q_i = Q_g^\mathrm{before}$$

The resulting square linear system is solved by Gaussian elimination with partial pivoting, and the conductor elements’ q_elem values are replaced.

11.4 Dielectric metadata

surface_model="dielectric" and epsilon_r are metadata in the current version. Dielectric polarization and dielectric boundary conditions are not yet reflected in the field calculation.

---

12. Statistics, history, and restart

Source: bem_simulator_stats, bem_simulator_io, bem_output_writer, bem_restart

12.1 Batch outcomes

count_batch_outcomes counts local-rank particles in five categories.

Index	Meaning
1	Total particles in the batch
2	Particles absorbed by the mesh
3	Particles counted as escaped
4	Particles that left through an open box boundary
5	Particles that survived to max_step

Among particles without absorbed_flag, those removed by an open boundary are escaped_boundary; those still alive at the end are survived_max_step.

In MPI runs, the batch-count array is allreduced so particles on non-root ranks are included in statistics.

12.2 History output

When history_stride > 0, BEACH writes charge_history.csv. The output condition is:

$$(stats.batches - 1) \bmod history\_stride = 0$$

so batch 1 is always written.

When output.write_potential_history=true, potential_history.csv is written at the same stride. For potential history, the field solver is refreshed with the current q_elem, and centroid potentials are computed and written.

12.3 Final output

When output.write_files=true, the root rank writes the main final outputs:

• summary.txt
• charges.csv
• mesh_potential.csv when enabled
• mesh_triangles.csv
• mesh_sources.csv

All ranks save checkpoint RNG state and macro residuals. In MPI runs, file names are rank-local.

12.4 Restart consistency

On restart, BEACH validates:

• checkpoint mesh_nelem matches the current mesh element count
• MPI world size matches the previous run
• statistics in summary.txt are finite and non-negative
• charges.csv element count and charge values are valid
• RNG state and macro residuals can be loaded

If a required checkpoint is missing, output.resume=true stops instead of falling back to a new run.

---

13. Parallelization and performance profiling

Source: bem_mpi, bem_performance_profile, bem_simulator_loop

13.1 OpenMP

Particle tracking is parallelized over particle indices with OpenMP.

• dq_thread(nelem, nth) accumulates collision charge thread-locally.
• The schedule is dynamic, 1 to reduce load imbalance from different particle lifetimes.
• Some OpenMP loops are also used in field-solver refresh and treecode node accumulation.

13.2 MPI

MPI parallelism splits particle generation and tracking by rank. Each rank holds the mesh and q_elem; at batch commit, dq is allreduced so every rank has the same charge state.

Main allreduces:

• sum of dq(nelem)
• sum of batch outcome counts

Only the root rank writes human-readable final CSV and history files. RNG state and macro residuals are saved per rank.

13.3 Performance profile

Set BEACH_PROFILE=1 to write times for major regions to performance_profile.csv. Main regions:

• load_or_init
• field_solver_init
• prepare_batch
• field_refresh
• particle_batch
• commit_charge
• mpi_reduce
• stats_update
• history_write
• write_results
• write_checkpoint

For MPI scaling evaluation, it is usually best to inspect the maximum time across ranks, rank_max_s.

---