Particle collision and boundary events

BEACH advances a particle only to the first mesh collision or box-boundary event on the current trajectory segment. This page is the reference for intersection search, event ordering, remainder advancement, and fail-closed statuses.

Accept only the first collision or crossing

Given current state $(\mathbf{x}_0,\mathbf{v}_0)$ and candidate $(\mathbf{x}_1,\mathbf{v}_1)$ from the Boris particle update, BEACH proceeds as follows:

1. Query the first mesh hit on $\mathbf{x}(t)=\mathbf{x}_0+t(\mathbf{x}_1-\mathbf{x}_0)$.
2. If the endpoint is inside the box, commit either the mesh hit or the endpoint.
3. If it is outside, compare the first box-face fraction against the mesh-hit $t$.
4. Absorb when the mesh is simultaneous or earlier; otherwise apply the box-face action.
5. If reflect, redistributed_reflect, or periodic keeps the particle alive, rebuild a candidate for the remaining time.

If mesh-hit and box-face fractions differ by at most $64\epsilon_\mathrm{mach}\max(1,|t|)$, the mesh is treated as first. Even if one particle step contains several events, BEACH always advances to the first event of the current segment before re-integrating the remainder.

Earliest event	Committed action
mesh	Return hit position and element index; absorb the particle
open face	Apply particle_boundary.ordinary_open_model at the event position
reflect face	Reverse normal velocity and advance the remainder from just inside the box
redistributed_reflect face	Reverse normal velocity, uniformly redistribute in-plane position, and advance the remainder
periodic face	Move just inside the opposite face, retain velocity, and advance the remainder

Intersect the trajectory segment with triangles

Reduce the triangle candidates

Mesh initialization builds an axis-aligned bounding box (AABB) for every triangle.

Element count	Candidate search
below 64	Linearly inspect every element AABB
64 and above	Use a uniform grid and 3D DDA to inspect only cells crossed by the trajectory segment

The uniform grid targets eight elements per cell and caps each axis at 128 cells. Triangle indices are stored in CSR form in every cell overlapped by the triangle AABB. These are fixed implementation values, not input parameters.

At query time, BEACH intersects the trajectory segment with the grid AABB and visits crossed cells with 3D DDA. A triangle may be registered in more than one cell, but only the smallest intersection parameter is retained. The DDA iteration bound is nx + ny + nz + 3. If cell indices, increments, or parameters do not make bounded finite progress, the result is not treated as “no hit”; the query returns collision_query_grid_stalled.

Select the first intersection

Each candidate triangle uses Möller–Trumbore intersection. Writing the triangle as

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

an intersection is accepted only when

$$0\le u\le1, \qquad 0\le v, \qquad u+v\le1, \qquad 0\le t\le1.$$

If the determinant magnitude is at most $64\epsilon_\mathrm{mach}$ times the product of segment length and the two edge lengths, the geometry is considered degenerate or nearly parallel. The determinant sign is not culled, so collisions are detected from both sides of a triangle. Triangle winding does not determine collision sidedness; it is used separately for quantities such as the field vacuum-side trace.

The smallest $t$ is selected among multiple triangle hits, with $\mathbf{h}=\mathbf{x}_0+t(\mathbf{x}_1-\mathbf{x}_0)$.

Map periodic-image collisions to the primary cell

With periodic2, the mesh still stores only base elements in the primary cell. BEACH enumerates only image shifts whose canonical mesh AABB can overlap the trajectory-segment AABB on the two periodic axes. For period $L$, segment range $[p_{\min},p_{\max}]$, and mesh range $[m_{\min},m_{\max}]$ on one axis,

$$n_{\min}=\left\lceil \frac{p_{\min}-m_{\max}-\mathrm{tol}}{L} \right\rceil, \qquad n_{\max}=\left\lfloor \frac{p_{\max}-m_{\min}+\mathrm{tol}}{L} \right\rfloor.$$

For each image, the trajectory segment is shifted by $-nL$ into the base-mesh frame and passed to the ordinary intersection query. A hit stores:

Value	Meaning
hit%pos	Physical coordinate on the periodic image that was hit
hit%pos_wrapped	Coordinate wrapped into the primary cell
hit%image_shift	Image indices on the two periodic axes
hit%elem_idx	Base-mesh element index

If candidate $t$ values agree within a relative tolerance of 1e-12, selection is deterministic by element index, first image index, then second image index. If the image count on one axis or its Cartesian product exceeds 4,096, the query does not skip work and report no hit; it returns collision_query_image_limit.

This image range is determined by which mesh images a particle trajectory segment can hit. It is unrelated to field_periodic_image_layers, which controls the field sum. See periodic2 electrostatics for field images.

Apply the action of each box face

Both faces on an axis listed in domain.periodic_axes are periodic. On nonperiodic faces, [particle_boundary] keys x_low, x_high, y_low, y_high, z_low, and z_high select open, reflect, or redistributed_reflect. Particle-boundary tables cannot specify periodic. When a segment reaches several faces at an edge or corner, faces within a machine-epsilon tolerance of the minimum fraction are combined into one mask.

Open, reflect, redistributed_reflect, and periodic faces

With particle_boundary.ordinary_open_model="escape", any open face in the mask removes the particle and increments escaped_boundary. Both reflection actions reverse the corresponding velocity components; periodic faces move the particle to the opposite side without changing velocity. A survivor’s event-axis coordinate is placed at an inward guard scaled to the box.

Ordinary reflect preserves the event position’s tangential components. redistributed_reflect applies the same velocity reflection but relocates the position. For a single-face event, it uniformly resamples both in-plane coordinates over the box span excluding the guards at its ends. If a simultaneous edge or corner mask contains redistributed_reflect, only axes outside the event mask are uniformly resampled. Event axes are not resampled and are placed at their respective inward guards. A corner therefore has no coordinate to resample, while an edge resamples only its remaining axis.

In-plane samples do not use a shared random stream. They are derived from sim.rng_seed and counters for batch, rank, particle, step, event, and axis. A fixed seed and MPI layout are reproducible independently of OpenMP scheduling; identical trajectories are not guaranteed after changing the rank decomposition.

[particles.species.boundary] overrides the same six faces for one species with inherit, open, reflect, or redistributed_reflect. inherit uses global [particle_boundary]. Faces selected by domain.periodic_axes are topology and cannot be overridden globally or per species. See Photoelectron emission and lifecycle for the closed-PE combination.

Reflection and periodic actions at a corner are applied from one face mask, making the result independent of axis traversal order.

Potential barrier

particle_boundary.ordinary_open_model="potential_barrier" compares outward normal kinetic energy at a single open face,

$$K_n=\frac{1}{2}m v_\mathrm{out}^2,$$

against

$$\Delta U=q\left(\phi_\infty-\phi_\mathrm{boundary}\right).$$

The particle reflects when $v_\mathrm{out}>0$, $\Delta U>0$, and $K_n<\Delta U$; otherwise it escapes. This is a reduced single-face model. A corner involving multiple simultaneous open faces is not generalized and returns particle_step_ambiguous_open_corner.

Reservoir inflow correction and closed PE are outside this page’s scope. See reservoir injection and particle escape and local return.

Advance the time remaining after a boundary crossing

The position at a box event is obtained from the candidate chord fraction, and every coordinate in the event mask is set exactly to its corresponding box face. The event velocity follows the chord tangent. Its speed is reconstructed from the discrete work of the predicted-midpoint electric field:

$$\lVert\mathbf{v}_\mathrm{event}\rVert^2=\lVert\mathbf{v}_0\rVert^2+ 2(q/m)\mathbf{E}_\mathrm{mid}\cdot(\mathbf{x}_\mathrm{event}-\mathbf{x}_0).$$

Thus an outward chord crossing cannot be processed with an inward partial-step Boris velocity, and a pure magnetic field preserves the event-speed norm. A non-positive or non-finite reconstructed speed, or a zero-length chord, fails closed as an unresolved event. BEACH also uses the chord crossing fraction as the remaining-time fraction, but this is not the exact physical crossing time under acceleration or magnetic rotation. Check the effect of this time-discretization approximation on boundary actions and potential-barrier decisions with dt, dt/2, and, when needed, dt/4. After applying the surviving action, BEACH uses a guard derived from the box coordinates and span to place reflected or periodic event-axis coordinates inside the face. This avoids a subnormal one-ULP offset at a zero-valued face and prevents the next event fraction from underflowing to zero.

$$\Delta t_\mathrm{remain}=(1-t_\mathrm{event})\Delta t_\mathrm{segment}$$

to build a new predicted-midpoint field and Boris candidate. The earliest mesh or box event is then queried on the remainder. The pre-reflection field sample is not reused for that remaining time.

One local continuation processes at most eight box events. If a ninth is required, particle_step_multiple_box_events is returned and the incomplete state is not committed. The limit detects an excessively large dt, narrow box, or fast particle; changing the scale-aware guard does not change it.

With multiple_box_events_retry_backend="upper_panel_fourier", a cached_kneq0 configuration replays only this failed step from its original position and velocity. The retry expands triangle-P0 charge through periodic2.reference_mode_layers and factorizes the nonzero Fourier field into its exponentially decaying form above the maximum z coordinate of every mesh vertex. Each evaluation adds the same periodic zero mode and sim.e0 exactly once. If any retry field sample is outside this upper-vacuum domain, the replay fails. A potential-barrier event evaluates its boundary potential with the same expansion and potential gauge, and its sample must satisfy the same domain. If that sample is outside the domain, or if the replay still cannot complete the events, BEACH retains the original multiple_box_events status and applies the configured policy. This is not an outer-plasma or sheath model, and it does not replace the FMM backend for ordinary particles. The geometry response is built once when the snapshot is initialized; multiplication by the current panel charges occurs only during retry field and potential evaluations. Mode count and geometry-response memory scale quadratically with reference_mode_layers, so choose it by checking seam-error and replay-result convergence.

The default multiple_box_events_policy="abort" fails the run closed at this point. Only an explicit "soft_discard" removes the affected macro-particle. Standard error records only the per-batch, all-rank count and absolute macro charge. Let $D$ be the cumulative discard count, $P$ the cumulative number of macro-particles processed in accepted batches, and $Q$ the cumulative absolute macro charge. Before commit, BEACH stops when

$$\left(D>G\ \text{and}\ \frac{D}{P}>f_{\mathrm{limit}}\right).$$

$G$ is multiple_box_events_soft_discard_count_grace, not a standalone count limit, and $f_{\mathrm{limit}}$ is multiple_box_events_soft_discard_fraction_limit. Equality is allowed in each comparison. multiple_box_events_soft_discard_abs_charge_limit is a warning threshold, not a stop condition; BEACH reports the first crossing of the cumulative absolute charge. summary.txt and checkpoints retain multiple_box_events_soft_discarded, multiple_box_events_soft_discarded_abs_charge_C, and multiple_box_events_soft_discard_fraction derived from $D/P$; the charge ledger also records discarded charge. Because a long normal history can dilute a late burst in the cumulative fraction, also audit the per-batch aggregate log. This is a bounded numerical workaround for qualitative comparisons, not a replacement for a physical boundary model. Regardless of the fallback policy, summary.txt records replay counts as multiple_box_events_retry_attempted and multiple_box_events_retry_resolved.

Stop when the query cannot be completed

A collision query is ok only when all required candidates were examined.

Status	Code	Meaning
collision_query_ok	0	Query completed
collision_query_image_limit	1	Periodic image enumeration exceeded 4,096
collision_query_index_range	2	Image bounds were non-finite or outside integer range
collision_query_invalid_segment	3	A trajectory-segment endpoint was non-finite
collision_query_grid_stalled	4	Invalid grid geometry or DDA failed to progress
particle_step_invalid_boundary	1001	Invalid particle, box, or event geometry
particle_step_multiple_box_events	1002	A ninth box event was needed in one step
particle_step_ambiguous_open_corner	1003	Multiple open faces occurred at a potential-barrier event

The former name particle_step_unsupported_barrier_corner remains as a compatibility alias for code 1003.

Treating these states as “no hit” could let particles pass through a surface, so tracking fails closed. OpenMP selects the smallest particle/step failure. MPI shares the failing rank and state so every rank reports the same batch/rank/particle/step/status. Photo raycasts apply the same rule to species/ray/bounce.

Set BEACH_COLLISION_DIAGNOSTICS=1 to inspect grid_stalled. BEACH writes the failing DDA branch, p0 / p1, grid bounds, cell indices, and values such as t_cur, t_next, and t_delta to standard error. The variable changes diagnostics only.

Converge collision positions and charging results

1. Halve dt and verify stability of the segment-hit element and position.
2. Build reduced cases with thin surfaces, edge and corner hits, and nearly parallel trajectory segments.
3. Test a mesh hit during the remainder after reflection and after periodic wrapping.
4. At a periodic seam, verify that pos and pos_wrapped identify the intended base element.
5. If multiple_box_events occurs, reduce dt first. With upper_panel_fourier, inspect its resolution rate and validity domain. With soft discard, verify that its rate and absolute charge cannot affect the conclusion.
6. Refine the mesh and check convergence of the first hit and final surface-charge distribution.

Code reference

• Box/mesh event ordering and remainder: bem_particle_stepper.f90
• Triangles, grid DDA, and periodic images: bem_collision.f90
• Box-face detection and actions: bem_boundary.f90
• Collision-grid construction: bem_mesh.f90
• Event ordering and eight-event-bound tests: test_particle_stepper.f90
• Boundary unit tests: test_boundary.f90
• Collision, periodic, and fail-closed tests: test_dynamics_basic.f90