simulation.cpp

#include <cmath>
#include <iostream>
 
#include <ostream>
#include <random>
#include <toml++/toml.hpp>
#include <Eigen/SparseCholesky>
 
#include "simulation.hpp"
#include "Eigen/Core"
#include "constants.hpp"
#include "interface.hpp"
#include "particles.hpp"
#include "poisson.hpp"
 
 
Simulation::Simulation(const IOInterface *interface) : 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);
 
};
 
void Simulation::run(){
 
    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;
};
 
 
void Simulation::postprocess(){
 
    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);
 
};