Main driver class coordinating IO, Poisson solve, particle push, and data storage for a 1D PIC simulation. More...

#include <simulation.hpp>

Public Member Functions
	Simulation (const IOInterface *interface)
	Construct the simulation from an IO interface.

void	run ()
	Execute the complete PIC time integration loop.

void	postprocess ()
	Write simulation results to output files.

Detailed Description

Main driver class coordinating IO, Poisson solve, particle push, and data storage for a 1D PIC simulation.

Definition at line 16 of file simulation.hpp.

Constructor & Destructor Documentation

◆ Simulation()

Simulation ( const IOInterface * interface )

explicit

Construct the simulation from an IO interface.

Parameters

interface Pointer to configuration and IO handler.

Definition at line 17 of file simulation.cpp.

                                                   : m_interface(interface) {
 
    std::cout << " -> Initializing simulation ..." << std::endl;
    
    // reproducable "random" numbers
    m_rng.seed(12345);
 
    // Init Poisson field solver
    Eigen::Vector4d pos = m_interface->load_vector<Eigen::Vector4d>("geometry.positions");
    Eigen::Vector4d pots = m_interface->load_vector<Eigen::Vector4d>("geometry.potentials");
    Eigen::Vector3i cells = m_interface->load_vector<Eigen::Vector3i>("discretization.space.cells");
    m_poisson = new PoissonSolver(pos, pots, cells);
 
    // Init particle storage
    double sp_weight = m_interface->load_scalar<double>("solve.particles.weight");
    m_particles = new Particles(0, constants::q_e, constants::m_e, sp_weight);
 
    // Init solution storage
    double t_final = m_interface->load_scalar<double>("discretization.time.t_final");
    double dt = m_interface->load_scalar<double>("discretization.time.dt");
    m_cellnumber = cells.sum();
    m_timesteps = static_cast<int>(std::ceil(t_final / dt));    
 
    m_time = Eigen::VectorXd::LinSpaced(m_timesteps, 0, m_timesteps * dt);
    m_phi = Eigen::MatrixXd::Zero(m_timesteps, m_cellnumber+1);
    m_rho = Eigen::MatrixXd::Zero(m_timesteps, m_cellnumber+1);
    m_efield = Eigen::MatrixXd::Zero(m_timesteps, m_cellnumber+1);
    
    // Init postproc storage
    m_bins = m_interface->load_scalar<int>("discretization.energy.bins");
    m_energy_dist = Eigen::MatrixXd::Zero(m_timesteps, m_bins);
    m_fluxes = Eigen::MatrixX4d::Zero(m_timesteps, 4);
 
};

Member Function Documentation

◆ postprocess()

void postprocess ( )

Write simulation results to output files.

Definition at line 116 of file simulation.cpp.

                            {
 
    std::cout << " -> Postprocessing results ..." << std::endl;
 
    // export results
    std::string comment = " | Time discretization";
    m_interface->write2csv("time_vec.csv", m_time, comment);
 
    comment = " | Cell sizes";
    m_interface->write2csv("dx_vec.csv", m_poisson->m_dx, comment);
    comment = " | Node positions";
    m_interface->write2csv("nodes_vec.csv", m_poisson->m_nodes, comment);
 
    comment = " | Charge density on nodes (rows: timesteps, cols: nodes)";
    m_interface->write2csv("rho_mat.csv", m_rho, comment);
    comment = " | Electric potential on nodes (rows: timesteps, cols: nodes)";
    m_interface->write2csv("phi_mat.csv", m_phi, comment);
    comment = " | Electric field in cells (rows: timesteps, cols: cells)";
    m_interface->write2csv("efield_mat.csv", m_efield, comment);
    comment = " | Electrode currents (rows: timesteps, cols: {1-Fil, 2-CE, 3-AG, 4-IC})";
    m_interface->write2csv("currents_mat.csv", m_fluxes, comment);
    comment = " | Electron kinetic energy distribution (rows: timesteps, cols: equidist. bins)";
    m_interface->write2csv("energy_mat.csv", m_energy_dist, comment);
 
};

◆ run()

void run ( )

Execute the complete PIC time integration loop.

Definition at line 52 of file simulation.cpp.

                    {
 
    std::cout << " -> Running the simulation ..." << std::endl;
    // general
    
    double source_flux_avg = m_interface->load_scalar<double>("solve.particles.source_flux");
    double initial_energy = m_interface->load_scalar<double>("solve.particles.init_energy");
    // collsions
    double p_ce = m_interface->load_scalar<double>("solve.collisions.prob_ce");
    double p_ag = m_interface->load_scalar<double>("solve.collisions.prob_ag");
    // max histogram energy
    double energy_max = m_interface->load_scalar<double>("discretization.energy.max");
 
    
    // Electron source with poisson distribution
    double dt = m_time[1] - m_time[0];
    std::poisson_distribution<int> pois(source_flux_avg * dt / m_particles->m_weight);
    double N_emit; // average number of electrons to emit over each timestep
 
    // Initial laplace solution
    std::cout << "[0] Solving initial Laplace " << std::endl;
    Eigen::VectorXd rho = Eigen::VectorXd::Zero(m_cellnumber+1);
    Eigen::VectorXd phi = m_poisson->solve_potential(rho);
    Eigen::VectorXd efield = m_poisson->solve_efield(phi);
    Eigen::VectorXd energies = Eigen::VectorXd::Zero(m_bins);
    Eigen::Vector4d fluxes = Eigen::Vector4d::Zero(4);
    m_phi.row(0) = phi.transpose();
    m_efield.row(0) = efield.transpose();
 
    // loop over rest of timestep here ...
    for (int i = 1; i < m_timesteps; i++) {
 
        std::cout << "["<< i << "] Solving timestep " << std::endl;
 
        // Emit on x=0 with boltzmann distributed velocities and poisson dist particle count
        N_emit = pois(m_rng);
        m_particles->emit(
            N_emit, 0, dt, initial_energy,
            m_poisson->m_nodes, efield, m_rng
            );
        // Move particles in efield solution
        m_particles->push_leapfrog(dt, m_poisson->m_nodes, efield);
        // Check for absorbtion
        fluxes = m_particles->fluxes(m_poisson->m_cap_xpos, m_rng, p_ce, p_ag);
        energies = m_particles->energy_distribution(m_bins, energy_max);
        // Deposite charge on grid
        rho = m_particles->deposit_charge(m_poisson->m_nodes, m_poisson->m_dx);
        // Update field solution w. Poisson solver
        phi = m_poisson->solve_potential(rho);
        efield = m_poisson->solve_efield(phi);
 
        // Save step quantities 
        m_rho.row(i) = rho.transpose();
        m_phi.row(i) = phi.transpose();
        m_efield.row(i) = efield.transpose();
        m_fluxes.row(i) = fluxes.transpose();
        m_energy_dist.row(i) = energies.transpose();
        
    }
    
    std::cout << " Finished simulation run " << std::endl;
};

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