poisson.cpp

#include <iostream>
#include <Eigen/SparseCholesky>
 
#include "poisson.hpp"
#include "Eigen/Core"
#include "constants.hpp"
 
PoissonSolver::PoissonSolver(Eigen::Vector4d positions, Eigen::Vector4d potentials, Eigen::Vector3i cells) : 
    m_cap_xpos(positions), m_cap_pots(potentials), m_cap_cells(cells) {
 
    // Discretizations (nodes and cell distances)
    m_nodes.resize(m_cap_cells.sum() + 1);
    m_nodes << 
        Eigen::VectorXd::Constant(1, m_cap_xpos[0]),
        Eigen::VectorXd::LinSpaced(m_cap_cells[0] + 1,m_cap_xpos[0],m_cap_xpos[1]).tail(m_cap_cells[0]),
        Eigen::VectorXd::LinSpaced(m_cap_cells[1] + 1,m_cap_xpos[1],m_cap_xpos[2]).tail(m_cap_cells[1]),
        Eigen::VectorXd::LinSpaced(m_cap_cells[2] + 1,m_cap_xpos[2],m_cap_xpos[3]).tail(m_cap_cells[2]);
    m_dx.resize(m_nodes.size() - 1);
    for (int i = 0; i < m_nodes.size() - 1; i++) {
        m_dx(i) = m_nodes(i+1) - m_nodes(i);
    }
    std::cout << "    Space discretization:" << 
        m_dx[0] << " | " << 
        m_dx[m_cap_cells[0]+1] << " | " << 
        m_dx[m_cap_cells[0]+m_cap_cells[1]+1] << std::endl; 
 
    // build finite difference matrices for each capacitor
    m_cap0_mat = this->tridiag(m_cap_cells[0] - 1);
    m_cap1_mat = this->tridiag(m_cap_cells[1] - 1);
    m_cap2_mat = this->tridiag(m_cap_cells[2] - 1);
 
};
 
Eigen::SparseMatrix<double> PoissonSolver::tridiag(int size){
 
    // sparse system matrix for second derivative in space
    typedef Eigen::Triplet<double> T;
    std::vector<T> triplets;
 
    for (int i = 0; i < size; ++i) {
        triplets.emplace_back(i, i, 2);      // main diagonal
        if (i > 0) triplets.emplace_back(i, i - 1, -1); // subdiagonal
        if (i < size-1) triplets.emplace_back(i, i + 1, -1); // superdiagonal
    }
 
    Eigen::SparseMatrix<double> A(size, size);
    A.setFromTriplets(triplets.begin(), triplets.end());
    return A;
 
};
 
 
Eigen::VectorXd PoissonSolver::solve_potential(Eigen::VectorXd& rho) {
 
    std::cout << " + solving poisson potential " << std::endl;
    double dx0 = m_dx[0];
    double dx1 = m_dx[m_cap_cells[0]+1];
    double dx2 = m_dx[m_cap_cells[0]+m_cap_cells[1]+1];
 
    // initial stage of the system is equal to laplace equation
    Eigen::VectorXd timestep_rhs0 = 
        rho.segment(1, m_cap_cells[0]-1) * std::pow(dx0,2) / constants::eps0;
    Eigen::VectorXd timestep_rhs1 = 
        rho.segment(m_cap_cells[0]+2, m_cap_cells[1]-1) * std::pow(dx1, 2) / constants::eps0;
    Eigen::VectorXd timestep_rhs2 = 
        rho.segment(m_cap_cells[0] + m_cap_cells[1] + 2, m_cap_cells[2]-1) * std::pow(dx2, 2) / constants::eps0;
    timestep_rhs0(0) += m_cap_pots[0];
    timestep_rhs0.tail(1)(0) += m_cap_pots[1];
    timestep_rhs1(0) += m_cap_pots[1];
    timestep_rhs1.tail(1)(0) += m_cap_pots[2];
    timestep_rhs2(0) += m_cap_pots[2];
    timestep_rhs2.tail(1)(0) += m_cap_pots[3];
 
    // solve system(s)
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver;
 
    solver.compute(m_cap0_mat);
    Eigen::VectorXd phi0 = solver.solve(timestep_rhs0);
    solver.compute(m_cap1_mat);
    Eigen::VectorXd phi1 = solver.solve(timestep_rhs1);
    solver.compute(m_cap2_mat);
    Eigen::VectorXd phi2 = solver.solve(timestep_rhs2);
    if(solver.info()!=Eigen::Success) {
        std::cerr << "Solving failed!" << std::endl;
    }
 
    // Merge potential solutions and include electrodes
    Eigen::VectorXd nodal_phi(m_nodes.size());
    nodal_phi <<
        Eigen::VectorXd::Constant(1, m_cap_pots[0]),
        phi0,
        Eigen::VectorXd::Constant(1, m_cap_pots[1]),
        phi1,
        Eigen::VectorXd::Constant(1, m_cap_pots[2]),
        phi2,
        Eigen::VectorXd::Constant(1, m_cap_pots[3]);
    
    return nodal_phi;
};
 
Eigen::VectorXd PoissonSolver::solve_efield(Eigen::VectorXd& nodal_phi){
    
    // 0th order e field calculation
    int n_nodes = nodal_phi.size();
    Eigen::VectorXd nodal_efield = Eigen::VectorXd::Zero(n_nodes);
 
    double hL;
    double hR;
    // ---------------------------
    // Central differences interior nodes
    // ---------------------------
    for (int i = 1; i < n_nodes-1; i++) {
        hL = m_dx[i-1];
        hR = m_dx[i];
 
        nodal_efield[i] = 
            (hL*hL * (nodal_phi[i+1]-nodal_phi[i]) 
            + hR*hR * (nodal_phi[i]-nodal_phi[i-1]))
            / (hL*hR * (hL+hR));
    }
 
 
    // left boundary (i = 0):
    {
        double h1 = m_dx[0];
        double h2 = m_dx[1];
        nodal_efield[0] = -(
              - (2*h1 + h2)/(h1*(h1 + h2)) * nodal_phi[0]
              + (h1 + h2)/(h1*h2) * nodal_phi[1]
              - h1/(h2*(h1 + h2)) * nodal_phi[2]
        );
    }
    // right boundary (i = N-1):
    {
        double h1 = m_dx[n_nodes-2];
        double h2 = m_dx[n_nodes-3];
        nodal_efield[n_nodes-1] = -(
              + (2*h1 + h2)/(h1*(h1 + h2)) * nodal_phi[n_nodes-1]
              - (h1 + h2)/(h1*h2) * nodal_phi[n_nodes-2]
              + h1/(h2*(h1 + h2)) * nodal_phi[n_nodes-3]
        );
    }
 
    return nodal_efield;
};