particles.cpp

#include "particles.hpp"
#include "constants.hpp"
 
#include <Eigen/Core>
#include <cmath>
#include <iostream>
#include <random>
 
// Constructors
// -----------------------------------------------------------------------------
Particles::Particles(int N, double q_in, double m_in, double w_in)
    : m_positions(N), m_vel_half(N), m_pos_prev(N), 
      m_charge(q_in), m_mass(m_in), m_weight(w_in) { }
 
// -----------------------------------------------------------------------------
// Resize positions + velocities
// -----------------------------------------------------------------------------
void Particles::resize(int N) {
    m_pos_prev.resize(N);
    m_positions.resize(N);
    m_vel_half.resize(N);
}
 
int Particles::find_cell_index(const Eigen::VectorXd& node_pos, double x) const {
    // Check if x is inside the domain
    if (x < node_pos[0] || x > node_pos[node_pos.size() - 1]) {
        throw std::out_of_range("Particle position " + std::to_string(x) + 
                                " is outside the domain [" +
                                std::to_string(node_pos[0]) + ", " +
                                std::to_string(node_pos[node_pos.size()-1]) + "]");
    }
 
    // Find the index of the cell containing x
    // Assumes node_pos is sorted ascending
    for (int i = 0; i < node_pos.size() - 1; ++i) {
        if (x >= node_pos[i] && x < node_pos[i + 1]) {
            return i;
        }
    }
 
    // Special case: x == last node
    return node_pos.size() - 2;
}
 
double Particles::interp_linear(const Eigen::VectorXd& x_nodes, const Eigen::VectorXd& val_nodes, double x) {
    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);
}
 
 
// -----------------------------------------------------------------------------
// Maxwellian velocity initialization
// T in Joules (k_B*T) or use T[eV] * e
// -----------------------------------------------------------------------------
void Particles::emit(
    int N, double xpos, double dt, double energy_eV, 
    const Eigen::VectorXd& node_pos, const Eigen::VectorXd& efield, std::mt19937& rng) {
 
    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;
 
            }
        }
    }
}
 
// -----------------------------------------------------------------------------
// Particle push (simple leapfrog / Euler method)
// E_part has the electric field interpolated TO particle positions
// -----------------------------------------------------------------------------
void Particles::push_leapfrog(
    double dt, const Eigen::VectorXd& node_pos, const Eigen::VectorXd& efield) {
    
    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;
    }
 
}
 
Eigen::Vector4d Particles::fluxes(const Eigen::Vector4d& el_pos, std::mt19937& rng, double p_ce, double p_ag) {
 
    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;
}
 
Eigen::VectorXd Particles::deposit_charge(const Eigen::VectorXd& node_pos,
                                          const Eigen::VectorXd& dx) const {
    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;
}
 
Eigen::VectorXd Particles::get_kinetic_energies() const {
    // kinetic E = 0.5 * m * v^2 / e   (eV) for physical particle
    Eigen::VectorXd energies = 0.5 * m_mass * (m_vel_half.array().square()) / constants::q_e;
    std::cout << " max v_p " << m_vel_half.maxCoeff() / constants::c << " (fraction of c0)" << std::endl;
    std::cout << " max E_kin " << energies.maxCoeff() << " (eV)" << std::endl;
    return energies;
}
 
Eigen::VectorXd Particles::energy_distribution(int n_bins, double W_max) const {
    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;
}