Particles Class Reference

Stores and advances a dynamically sized set of 1D superparticles. More...

#include <particles.hpp>

Public Member Functions
	Particles ()=default
	Default constructor, creates an empty particle set.
	Particles (int N, double q, double m, double w)
	Construct a particle set with N superparticles.
void	resize (int N)
	Resize particle arrays to contain N superparticles.
double	interp_linear (const Eigen::VectorXd &x_nodes, const Eigen::VectorXd &E_nodes, double x)
void	emit (int N, double xpos, double dt, double energy_eV, const Eigen::VectorXd &node_pos, const Eigen::VectorXd &efield, std::mt19937 &rng)
	Emit N particles at position xpos with Maxwellian velocity distribution.
void	push_leapfrog (double dt, const Eigen::VectorXd &node_pos, const Eigen::VectorXd &efield)
	Advance particles using leapfrog (velocity half-step scheme).
Eigen::Vector4d	fluxes (const Eigen::Vector4d &el_pos, std::mt19937 &rng, double p_ce, double p_ag)
	Compute absorbed particle fluxes on four electrodes.
Eigen::VectorXd	deposit_charge (const Eigen::VectorXd &node_pos, const Eigen::VectorXd &dx) const
	Deposit particle charge onto grid nodes (CIC scheme).
Eigen::VectorXd	energy_distribution (int n_bins, double W_max) const

Data Fields
Eigen::VectorXd	m_positions
	Particle positions (1D).
Eigen::VectorXd	m_pos_prev
	Previous positions for boundary‐crossing detection.
Eigen::VectorXd	m_vel_half
	Half-step velocities for leapfrog integration.
double	m_charge
	Charge of a single physical particle (C).
double	m_mass
	Mass of a single physical particle (kg).
double	m_weight = 1.0
	Superparticle statistical weight.

Detailed Description

Stores and advances a dynamically sized set of 1D superparticles.

Handles particle emission, leapfrog pushing under electric fields, charge deposition onto grid nodes, and absorbing boundary flux diagnostics.

Definition at line 13 of file particles.hpp.

Constructor & Destructor Documentation

Particles()

Particles	(	int	N,
		double	q,
		double	m,
		double	w )

Construct a particle set with N superparticles.

Parameters

N	Number of particles.
q	Charge per physical particle.
m	Mass per physical particle.
w	Statistical weight per superparticle.

Definition at line 11 of file particles.cpp.

    : m_positions(N), m_vel_half(N), m_pos_prev(N), 
      m_charge(q_in), m_mass(m_in), m_weight(w_in) { }

Member Function Documentation

deposit_charge()

Eigen::VectorXd deposit_charge	(	const Eigen::VectorXd &	node_pos,
		const Eigen::VectorXd &	dx ) const

Deposit particle charge onto grid nodes (CIC scheme).

Parameters

node_pos	Node positions.
dx	Cell widths.

Returns

Charge density vector (C/m) at nodes.

Definition at line 216 of file particles.cpp.

                                                                         {
    std::cout << " + depositing particle charges (OMP)" << std::endl;
 
    int Nnodes = node_pos.size();
    int Np     = m_positions.size();
 
    Eigen::VectorXd rho = Eigen::VectorXd::Zero(Nnodes);
 
    double qsp = m_charge * m_weight;
 
    #pragma omp parallel
    {
        // Thread-local accumulator
        Eigen::VectorXd local_rho = Eigen::VectorXd::Zero(Nnodes);
 
        #pragma omp for nowait
        for (int p = 0; p < Np; p++) {
 
            double xpos = m_positions[p];
            int j = find_cell_index(node_pos, xpos);
 
            double xL = node_pos[j];
            double xR = node_pos[j+1];
 
            double wL = (xR - xpos) / dx[j];
            double wR = (xpos - xL) / dx[j];
 
            local_rho[j]     += qsp * wL;
            local_rho[j + 1] += qsp * wR;
        }
 
        // reduction (atomic add)
        #pragma omp critical
        {
            rho += local_rho;
        }
    }
 
    // Convert charge to density
    #pragma omp parallel for
    for (int j = 0; j < Nnodes - 1; j++) {
        if (dx[j] > 0) {
            rho[j]     /= dx[j];
            rho[j + 1] /= dx[j];
        }
    }
 
    return rho;
}

emit()

void emit	(	int	N,
		double	xpos,
		double	dt,
		double	energy_eV,
		const Eigen::VectorXd &	node_pos,
		const Eigen::VectorXd &	efield,
		std::mt19937 &	rng )

Emit N particles at position xpos with Maxwellian velocity distribution.

Parameters

N	Number of emitted particles.
xpos	Injection position.
dt	Injection time intervall.
energy_eV	Thermal energy (eV) for Maxwellian sampling.
rng	Random generator for velocity/time sampling.

Definition at line 71 of file particles.cpp.

                                                                                 {
 
    std::cout << " + Emmiting " << N << " new particles" << std::endl;
    int n_nodes = node_pos.size();
    int old_N = m_positions.size();
    int new_N = old_N + N;
 
    m_pos_prev.conservativeResize(new_N);
    m_positions.conservativeResize(new_N);
    m_vel_half.conservativeResize(new_N);
 
    double kT = energy_eV * constants::q_e;
    double sigma_v = std::sqrt(kT / m_mass);
 
    double t_emit, v_emit = 0;
    double a = (m_charge * efield[0]) / m_mass;
    bool overshoot = false;
 
    // IMPORTANT: RNG must be thread-local
    std::normal_distribution<double> normal(0.0, sigma_v);  // Maxwellian velocity
    std::uniform_real_distribution<double> uniform(0.0, dt);        // Uniform emission time
 
    #pragma omp parallel
    {
        std::mt19937 local_rng(rng());  // copy engine per thread
 
        #pragma omp for
        for (int i = 0; i < N; i++) {
            int idx = old_N + i;
 
            t_emit = uniform(local_rng); // relative time of emission inside timestep
            v_emit =  std::fabs(normal(local_rng)); // thermionic emission speed
            double p_pos = xpos; 
            double p_vel = v_emit;
 
            // initial values
            m_pos_prev[idx] = xpos;
            m_positions[idx] = xpos;
            m_vel_half[idx] = v_emit;
            
            // preadvance position and velocity in efield solution with substepping
            double dt_sub = dt / std::ceil(n_nodes/2);
            double local_a = 0;
            for (int j = 0; j < std::ceil(n_nodes/2); j++) {
                local_a = (m_charge * this->interp_linear(node_pos, efield, p_pos)) / m_mass;
                m_vel_half[idx] += local_a * dt_sub;
                m_positions[idx] += m_vel_half[idx] * dt_sub;
 
            }
        }
    }
}

energy_distribution()

Eigen::VectorXd energy_distribution	(	int	n_bins,
		double	W_max ) const

Definition at line 275 of file particles.cpp.

                                                                           {
    Eigen::VectorXd energies = get_kinetic_energies();
 
    Eigen::VectorXd hist = Eigen::VectorXd::Zero(n_bins);
    double dW = W_max / n_bins;
 
    for (int i = 0; i < energies.size(); ++i) {
        int bin = std::min(static_cast<int>(energies[i] / dW), n_bins - 1);
        hist(bin) += m_weight;  // scale by superparticle weight
    }
 
    // normalize if desired (number of particles per unit energy)
    // hist /= dW;
    return hist;
}

fluxes()

Eigen::Vector4d fluxes	(	const Eigen::Vector4d &	el_pos,
		std::mt19937 &	rng,
		double	p_ce,
		double	p_ag )

Compute absorbed particle fluxes on four electrodes.

Flux order: 0: filament
1: control electrode
2: acceleration grid
3: ion collector

Randomized absorption at the control electrode and acceleration grid is controlled by probabilities p_ce and p_ag.

Parameters

el_pos	Positions of the four electrodes.
rng	Random generator for absorption events.
p_ce	Absorption probability at control electrode.
p_ag	Absorption probability at acceleration grid.

Returns

Vector of absorbed charges per boundary.

Definition at line 153 of file particles.cpp.

                                                                                                      {
 
    int N = m_positions.size();
    int initial_size = N;
 
    int fil_particles = 0;
    int ce_particles  = 0;
    int ag_particles  = 0;
    int ic_particles  = 0;
 
    int i = 0;
    while (i < N) {
 
        bool removed = false;
 
        if (m_positions[i] < el_pos[0]) {
            // --- Particle hit filament ---
            fil_particles++;
            removed = true;
        } 
        else if (m_positions[i] > el_pos.tail(1).value()) {
            // --- Particle hit ion collector ---
            ic_particles++;
            removed = true;
        } 
        else if ((m_pos_prev[i] - el_pos[1]) * (m_positions[i] - el_pos[1]) <= 0.0) {
            if (std::generate_canonical<double, 53>(rng) < p_ce) {
                ce_particles++;
                removed = true;
            }
        }
        else if ((m_pos_prev[i] - el_pos[2]) * (m_positions[i] - el_pos[2]) <= 0.0) {
            // --- Particle crossed acceleration grid ---
            if (std::generate_canonical<double, 53>(rng) < p_ag) {
                ag_particles++;
                removed = true;
            }
        }
 
        if (removed) {
            // Swap-and-pop to remove particle efficiently
            int last = N - 1;
            m_positions[i] = m_positions[last];
            m_pos_prev[i] = m_pos_prev[last];
            m_vel_half[i] = m_vel_half[last];
            --N;  // reduce particle count
        } 
        else {
            ++i;
        }
    }
 
    // Resize vectors to remove absorbed particles
    m_positions.conservativeResize(N);
    m_pos_prev.conservativeResize(N);
    m_vel_half.conservativeResize(N);
 
    Eigen::Vector4d particle_fluxes(fil_particles, ce_particles, ag_particles, ic_particles);
 
    std::cout << " + removed " << initial_size - N << " particles " << std::endl;
    return particle_fluxes * m_charge * m_weight;
}

interp_linear()

double interp_linear	(	const Eigen::VectorXd &	x_nodes,
		const Eigen::VectorXd &	E_nodes,
		double	x )

Definition at line 45 of file particles.cpp.

                                                                                                      {
    const int N = x_nodes.size();
 
    // --- bounds check (clamp) ---
    if (x <= x_nodes[0])     return val_nodes[0];
    if (x >= x_nodes[N - 1]) return val_nodes[N - 1];
 
    // --- locate cell via binary search ---
    auto it = std::lower_bound(x_nodes.data(), x_nodes.data() + N, x);
    int i = std::max(0, int(it - x_nodes.data() - 1));
 
    // --- linear interpolation ---
    double x0 = x_nodes[i];
    double x1 = x_nodes[i+1];
    double E0 = val_nodes[i];
    double E1 = val_nodes[i+1];
 
    double t = (x - x0) / (x1 - x0);
    return E0 + t * (E1 - E0);
}

push_leapfrog()

void push_leapfrog	(	double	dt,
		const Eigen::VectorXd &	node_pos,
		const Eigen::VectorXd &	efield )

Advance particles using leapfrog (velocity half-step scheme).

Parameters

dt	Time step.
efield	Cell averaged electric fields.
node_pos	Nodal positions.

Definition at line 130 of file particles.cpp.

                                                                           {
    
    int Np = m_positions.size();
    m_pos_prev = m_positions;  // store previous positions for diagnostics
 
    #pragma omp parallel for
    for (int i = 0; i < Np; i++) {
        // 1. Interpolate electric field at particle position
        double local_efield = interp_linear(node_pos, efield, m_positions[i]);
 
        // 2. Compute acceleration
        double local_a = m_charge * local_efield / m_mass;
 
        // 3. Update half-step velocity
        m_vel_half[i] += local_a * dt;
 
        // 4. Update position
        m_positions[i] += m_vel_half[i] * dt;
    }
 
}

resize()

void resize

(

int

N

)

Resize particle arrays to contain N superparticles.

Parameters

N	New total particle count.

Definition at line 18 of file particles.cpp.

                            {
    m_pos_prev.resize(N);
    m_positions.resize(N);
    m_vel_half.resize(N);
}

Field Documentation

m_charge

double m_charge

Charge of a single physical particle (C).

Definition at line 35 of file particles.hpp.

m_mass

double m_mass

Mass of a single physical particle (kg).

Definition at line 38 of file particles.hpp.

m_pos_prev

Eigen::VectorXd m_pos_prev

Previous positions for boundary‐crossing detection.

Definition at line 29 of file particles.hpp.

m_positions

Eigen::VectorXd m_positions

Particle positions (1D).

Definition at line 26 of file particles.hpp.

m_vel_half

Eigen::VectorXd m_vel_half

Half-step velocities for leapfrog integration.

Definition at line 32 of file particles.hpp.

m_weight

double m_weight = 1.0

Superparticle statistical weight.

Definition at line 41 of file particles.hpp.

---

The documentation for this class was generated from the following files:

• src/particles.hpp
• src/particles.cpp