Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <cstdio>
- #include <cstdlib>
- #include <cstring>
- #include <cstdbool>
- #include <cmath>
- #include <ctime>
- #include <random>
- #include <vector>
- #include <string>
- #include <fstream>
- #include <sstream>
- #include <algorithm>
- #include <stdexcept>
- #include <iostream>
- #include <iomanip>
- using namespace::std;
- // constants
- const double PI = 3.141592653589793; // mathematical constant Pi
- const double TWO_PI = 2.0 * PI; // two times Pi
- const double E_CHARGE = 1.60217662e-19; // electron charge [C]
- const double EV_TO_J = E_CHARGE; // eV <-> Joule conversion factor
- const double E_MASS = 9.109e-31; // mass of electron [kg]
- const double HE_MASS = 6.67e-27; // mass of He atom [kg]
- const double MU_HEHE = HE_MASS / 2.0; // reduced mass of two He atoms [kg]
- const double K_BOLTZMANN = 1.38064852e-23; // Boltzmann's constant [J/K]
- const double EPSILON0 = 8.85418781e-12; // permittivity of free space [F/m]
- const double C_light = 299792458.0; // speed of light [m/s]
- // simulation parameters
- const int N_G = 501; // number of grid points
- const int N_T = 2000; // time steps within an RF period
- const double FREQUENCY = 13.56e6; // driving frequency [Hz]
- const double VOLTAGE = 350.0; // voltage amplitude [V]
- const double L = 0.04; // electrode gap [m]
- const double PRESSURE = 13.3322; // gas pressure [Pa] // n*k*T to match Turner's case
- const double T_neutral = 300.0; // background gas temperature [K] (also ion temperature)
- const double T_electron = 30000.0; // initial electron temperatutre [K]
- const double WEIGHT = 41875.0; // weight of superparticles
- const double ELECTRODE_AREA = 1.6e-4; // (fictive) electrode area [m^2]
- const int N_INIT = 65536; // number of initial electrons and ions
- // additional (derived) constants
- const double PERIOD = 1.0 / FREQUENCY; // RF period length [s]
- const double DT_E = PERIOD / (double)(N_T); // electron time step [s]
- const int N_SUB = 10; // ions move only in these cycles (subcycling)
- const int N_avg = 5; // number of cycles to average rates --- Artem
- int counter = 0; // cycles since update of rates --- Artem
- const double DT_I = N_SUB * DT_E; // ion time step [s]
- const double DX = L / (double)(N_G - 1); // spatial grid division [m]
- const double INV_DX = 1.0 / DX; // inverse of spatial grid size [1/m]
- const double GAS_DENSITY = PRESSURE / (K_BOLTZMANN * T_neutral); // background gas density [1/m^3]
- const double OMEGA = TWO_PI * FREQUENCY; // angular frequency [rad/s]
- const double Diff_coeff = 8.992E-6*pow(T_neutral,1.5)/(PRESSURE/133.3); // diffusion coefficient for metastable states
- // electron and ion cross sections
- const int N_CS = 20; // total number of processes / cross sections
- const int E_ELA = 0; // process identifier: electron/elastic
- const int E_EXC_1 = 1; // process identifier: electron/excitation1
- const int E_EXC_2 = 2; // process identifier: electron/excitation1
- const int E_ION = 3; // process identifier: electron/ionization
- const int I_ISO = 4; // process identifier: ion/elastic/isotropic
- const int I_BACK = 5; // process identifier: ion/elastic/backscattering
- const int E_SUPER_1 = 6; // 1^1S -> 2^3S
- const int E_SUPER_2 = 7; // 1^1S -> 2^1S
- const int E_EXC_3 = 8; // 1^1S -> 2^1P
- const int E_EXC_4 = 9; // 1^1S -> 2^3P
- const int E_EXC_12 = 10; // 2^1S -> 2^1P
- const int E_EXC_13 = 11; // 2^3S -> 2^3P
- const int E_EXC_17 = 12; // 2^3S -> 2^1S
- const int E_EXC_19 = 13; // 2^1S -> 2^3P
- const int E_SUPER_3 = 14; // 1^1S -> 2^1P
- const int E_SUPER_4 = 15; // 1^1S -> 2^3P
- const int E_SUPER_12 = 16; // 2^1S -> 2^1P
- const int E_SUPER_13 = 17; // 2^3S -> 2^3P
- const int E_SUPER_17 = 18; // 2^3S -> 2^1S
- const int E_SUPER_19 = 19; // 2^1S -> 2^3P
- const double E_EXC_TH_1 = 19.82; // 1^1S -> 2^3S threshold [eV]
- const double E_EXC_TH_2 = 20.61; // 1^1S -> 2^1S threshold [eV]
- const double E_EXC_TH_3 = 21.218; // 1^1S -> 2^1P threshold [eV]
- const double E_EXC_TH_4 = 20.96; // 1^1S -> 2^3P threshold [eV]
- const double E_EXC_TH_12 = 0.602; // 2^1S -> 2^1P threshold [eV]
- const double E_EXC_TH_13 = 1.1444; // 2^3S -> 2^3P threshold [eV]
- const double E_EXC_TH_17 = 0.796; // 2^3S -> 2^1S threshold [eV]
- const double E_EXC_TH_19 = 0.3483; // 2^1S -> 2^3P threshold [eV]
- const double E_ION_TH = 24.587; // electron impact ionization threshold [eV]
- const int CS_RANGES = 1000000; // number of entries in cross section arrays
- const double DE_CS = 0.001; // energy division in cross section arrays [eV]
- typedef float cross_section[CS_RANGES]; // cross section array
- cross_section sigma[N_CS]; // set of cross section arrays
- cross_section sigma_tot_e; // total macroscopic cross section of electrons
- cross_section sigma_tot_i; // total macroscopic cross section of ions
- double nu_avg; // average nu for electrons thorugh 1 cycle
- // particle coordinates
- const int MAX_N_P = 10000000; // maximum number of particles (electrons / ions)
- typedef double particle_vector[MAX_N_P]; // array for particle properties
- int N_e = 0; // number of electrons
- int N_i = 0; // number of ions
- int N_e1 = 0; // number of singlet excited states
- int N_e2 = 0; // number of triplet excited states
- particle_vector x_e, vx_e, vy_e, vz_e; // coordinates of electrons (one spatial, three velocity components)
- particle_vector x_i, vx_i, vy_i, vz_i; // coordinates of ions (one spatial, three velocity components)
- typedef double xvector[N_G]; // array for quantities defined at gird points
- xvector efield,pot; // electric field and potential
- xvector e_density,i_density; // electron and ion densities
- xvector cumul_e_density,cumul_i_density; // cumulative densities
- //excited states handling
- xvector e1_density = {0.0}; // 2^3S state
- xvector e2_density = {0.0}; // 2^1S state
- xvector e3_density = {0.0}; // 2^1P state
- xvector e4_density = {0.0}; // 2^3P state
- xvector sum_electron_density = {0.0}; xvector avg_electron_density = {0.0};
- xvector sum_rate1f = {0.0}; xvector sum_rate1b = {0.0}; xvector sum_rate2f = {0.0}; xvector sum_rate2b = {0.0};
- xvector sum_rate3f = {0.0}; xvector sum_rate3b = {0.0}; xvector sum_rate4f = {0.0}; xvector sum_rate4b = {0.0};
- xvector sum_rate12f = {0.0}; xvector sum_rate12b = {0.0}; xvector sum_rate13f = {0.0}; xvector sum_rate13b = {0.0};
- xvector sum_rate17f = {0.0}; xvector sum_rate17b = {0.0}; xvector sum_rate19f = {0.0}; xvector sum_rate19b = {0.0};
- xvector avg_rate1f = {0.0}; xvector avg_rate1b = {0.0}; xvector avg_rate2f = {0.0}; xvector avg_rate2b = {0.0};
- xvector avg_rate3f = {0.0}; xvector avg_rate3b = {0.0}; xvector avg_rate4f = {0.0}; xvector avg_rate4b = {0.0};
- xvector avg_rate12f = {0.0}; xvector avg_rate12b = {0.0}; xvector avg_rate13f = {0.0}; xvector avg_rate13b = {0.0};
- xvector avg_rate17f = {0.0}; xvector avg_rate17b = {0.0}; xvector avg_rate19f = {0.0}; xvector avg_rate19b = {0.0};
- typedef unsigned long long int Ullong; // compact name for 64 bit unsigned integer
- Ullong N_e_abs_pow = 0; // counter for electrons absorbed at the powered electrode
- Ullong N_e_abs_gnd = 0; // counter for electrons absorbed at the grounded electrode
- Ullong N_i_abs_pow = 0; // counter for ions absorbed at the powered electrode
- Ullong N_i_abs_gnd = 0; // counter for ions absorbed at the grounded electrode
- // electron energy probability function
- const int N_EEPF = 2000; // number of energy bins in Electron Energy Probability Function (EEPF)
- const double DE_EEPF = 0.05; // resolution of EEPF [eV]
- typedef double eepf_vector[N_EEPF]; // array for EEPF
- eepf_vector eepf = {0.0}; // time integrated EEPF in the center of the plasma
- // ion flux-energy distributions
- const int N_IFED = 200; // number of energy bins in Ion Flux-Energy Distributions (IFEDs)
- const double DE_IFED = 1.0; // resolution of IFEDs [eV]
- typedef int ifed_vector[N_IFED]; // array for IFEDs
- ifed_vector ifed_pow = {0}; // IFED at the powered electrode
- ifed_vector ifed_gnd = {0}; // IFED at the grounded electrode
- double mean_i_energy_pow; // mean ion energy at the powered electrode
- double mean_i_energy_gnd; // mean ion energy at the grounded electrode
- // spatio-temporal (XT) distributions
- const int N_BIN = 20; // number of time steps binned for the XT distributions
- const int N_XT = N_T / N_BIN; // number of spatial bins for the XT distributions
- typedef double xt_distr[N_G][N_XT]; // array for XT distributions (decimal numbers)
- xt_distr pot_xt = {0.0}; // XT distribution of the potential
- xt_distr efield_xt = {0.0}; // XT distribution of the electric field
- xt_distr ne_xt = {0.0}; // XT distribution of the electron density
- xt_distr ni_xt = {0.0}; // XT distribution of the ion density
- xt_distr ue_xt = {0.0}; // XT distribution of the mean electron velocity
- xt_distr ui_xt = {0.0}; // XT distribution of the mean ion velocity
- xt_distr je_xt = {0.0}; // XT distribution of the electron current density
- xt_distr ji_xt = {0.0}; // XT distribution of the ion current density
- xt_distr powere_xt = {0.0}; // XT distribution of the electron powering (power absorption) rate
- xt_distr poweri_xt = {0.0}; // XT distribution of the ion powering (power absorption) rate
- xt_distr meanee_xt = {0.0}; // XT distribution of the mean electron energy
- xt_distr meanei_xt = {0.0}; // XT distribution of the mean ion energy
- xt_distr counter_e_xt = {0.0}; // XT counter for electron properties
- xt_distr counter_i_xt = {0.0}; // XT counter for ion properties
- xt_distr ioniz_rate_xt = {0.0}; // XT distribution of the ionisation rate
- xt_distr e1_xt = {0.0}; // XT distribution of 2^3S excited states density - Artem
- xt_distr e2_xt = {0.0}; // XT distribution of 2^1S excited states density - Artem
- xt_distr e3_xt = {0.0}; // XT distribution of 2^1P excited states density - Artem
- xt_distr e4_xt = {0.0}; // XT distribution of 2^3P excited states density - Artem
- double mean_energy_accu_center = 0; // mean electron energy accumulator in the center of the gap
- Ullong mean_energy_counter_center = 0; // mean electron energy counter in the center of the gap
- Ullong N_e_coll = 0; // counter for electron collisions
- Ullong N_i_coll = 0; // counter for ion collisions
- double Time; // total simulated time (from the beginning of the simulation)
- int cycle,no_of_cycles,cycles_done; // current cycle and total cycles in the run, cycles completed
- int arg1; // used for reading command line arguments
- char st0[80]; // used for reading command line arguments
- FILE *datafile; // used for saving data
- bool measurement_mode; // flag that controls measurements and data saving
- const double gamma_i = 0.3; // yielding coefficient for ions bombarding the surface
- const double gamma_m = 0.3; // yielding coefficient for metastables bombarding the surface - not using now
- const double gamma_p = -10.0; // yielding coefficient for photons bombarding the surface - not using now
- const double gamma_k = 0.5; // recombination coefficient
- const double r = 0.7; // probabilty for the electron to reflect at the surface
- const double T_SEE = 300; // initial temperature of SEE electrons [K]
- const double Ry = 13.6057; // Rydberg constant, eV
- const double Bohr = 0.8797e-20; // atomic consant related to Bohr Radius, m^2
- //---------------------------------------------------------------------------//
- // C++ Mersenne Twister 19937 generator //
- // R01(MTgen) will genarate uniform distribution over [0,1) interval //
- // RMB(MTgen) will generate Maxwell-Boltzmann distribution (of gas atoms) //
- //---------------------------------------------------------------------------//
- std::random_device rd{};
- std::mt19937 MTgen(rd());
- std::uniform_real_distribution<> R01(0.0, 1.0);
- std::normal_distribution<> RMB_n(0.0,sqrt(K_BOLTZMANN * T_neutral / HE_MASS));
- std::normal_distribution<> RMB_e(0.0,sqrt(K_BOLTZMANN * T_electron / E_MASS));
- std::normal_distribution<> RMB_SEE(0.0,sqrt(K_BOLTZMANN * T_SEE / E_MASS));
- double series(double tau){
- int k = 20; // accuracy of sum calculation, number of series expansion parts
- double sum = 0.0;
- for (int i = 1; i < k+1; i++){
- double factor = std::tgamma((i+2) + 1);
- // for (int j = 0; j < i+2; j++){
- // factor *= i + 2 - j;
- // }
- sum += pow(-1.0, i) * pow(tau, i+1) / (i*factor*sqrt(i+2));
- }
- return sum;
- }
- double compute_escape_factor(double tau){
- if (tau <= 0) return 1; //optically thin plasma
- if (tau > 2.5) return 1/(2.0*tau*sqrt(M_PI))*(sqrt(log(2.0*tau)) + 0.25/sqrt(log(2.0*tau)) + 0.14);
- else return 1.0 - 0.8293*tau + 0.7071*tau*log(2.0*tau) + series(tau);
- }
- class CSInterpolator {
- public:
- // load "filename" which must have two whitespace‐separated columns:
- // energy (eV) cross_section (in m^2 or cm^2 as you prefer)
- CSInterpolator(const std::string &filename) {
- std::ifstream in(filename);
- if (!in) throw std::runtime_error("CSInterpolator: cannot open " + filename);
- double E, sigma;
- std::string line;
- while (std::getline(in, line)) {
- std::istringstream iss(line);
- if (iss >> E >> sigma) {
- E_pts_.push_back(E);
- sigma_pts_.push_back(sigma);
- }
- }
- if (E_pts_.size()<2)
- throw std::runtime_error("CSInterpolator: need at least two data points in " + filename);
- }
- // return σ(E) by simple linear interpolation (clamped to end‐points)
- double operator()(double E) const {
- auto it = std::lower_bound(E_pts_.begin(), E_pts_.end(), E);
- if (it == E_pts_.begin()) {
- // std::cerr << "Warning: E="<<E<<" below data range, clamping to "<< 0.0 <<"\n";
- return 0.0;
- }
- if (it == E_pts_.end()){
- // std::cerr << "Warning: E="<<E<<" above data range, clamping to "<< sigma_pts_.back() <<"\n";
- return sigma_pts_.back();
- }
- size_t idx = (it - E_pts_.begin());
- double E1 = E_pts_[idx-1], E2 = E_pts_[idx];
- double s1 = sigma_pts_[idx-1], s2 = sigma_pts_[idx];
- // linear interp
- return s1 + (s2 - s1) * (E - E1)/(E2 - E1);
- }
- private:
- std::vector<double> E_pts_, sigma_pts_;
- };
- // Cross-Section function for excitation cross-sections.
- // here for type 'd' - dipole-forbidden, 's' - spin-forbidden, 'a' - dipole-allowed
- // g - stat weight of initial level (g = 1 for ground state). See paper: RalchenkoJanev2008
- double forbidden_CS(double A1, double A2, double A3, double A4, double A5, double A6, double E, double E_th, int g, char type){
- double x = E/E_th;
- double strength, S1, S2, CS;
- if (type == 'd'){
- S1 = (A1 + A2/x + A3/(x*x) + A4/(x*x*x));
- S2 = (x*x/(x*x + A5));
- }
- if (type == 's'){
- S1 = (A1 + A2/x + A3/(x*x) + A4/(x*x*x));
- S2 = (1.0/(x*x + A5));
- }
- if (type == 'a'){
- S1 = (A1*log(x) + A2 + A3/x + A4/(x*x) + A5/(x*x*x));
- S2 = (x + 1.0)/(x + A6);
- }
- strength = S1*S2;
- CS = Bohr * Ry/(g*E) * strength;
- return CS;
- }
- void set_electron_cross_sections_ar(void){
- int i;
- double A1_1, A1_2, A1_3, A1_4, A1_5, A1_6, A2_1, A2_2, A2_3, A2_4, A2_5, A2_6;
- double A3_1, A3_2, A3_3, A3_4, A3_5, A3_6, A4_1, A4_2, A4_3, A4_4, A4_5, A4_6;
- double A12_1, A12_2, A12_3, A12_4, A12_5, A12_6, A13_1, A13_2, A13_3, A13_4, A13_5, A13_6;
- double A17_1, A17_2, A17_3, A17_4, A17_5, A17_6, A19_1, A19_2, A19_3, A19_4, A19_5, A19_6;
- int g1, g2, g3, g4, g12, g13, g17, g19; // stat weights of initial level
- A1_1 = +6.888e-01; A1_2 = +1.975e-01; A1_3 = +7.232e+00; A1_4 = -4.839e+00; A1_5 = +5.003e+01; A1_6 = +0.000e+00; // 1¹S -> 2^3S -- s-f
- A2_1 = +1.888e-01; A2_2 = -5.754e-01; A2_3 = +3.439e+00; A2_4 = -2.088e+00; A2_5 = +2.544e+01; A2_6 = +0.000e+00; // 1¹S -> 2^1S -- d-f
- A3_1 = +7.087e-01; A3_2 = -9.347e-02; A3_3 = -1.598e+00; A3_4 = +2.986e+00; A3_5 = -1.293e+00; A3_6 = +3.086e-01; // 1¹S -> 2¹P -- d-a
- A4_1 = +2.823e-01; A4_2 = +2.048e+00; A4_3 = +5.287e+00; A4_4 = -7.363e+00; A4_5 = +2.728e+01; A4_6 = +0.000e+00; // 1¹S -> 2³P -- s-f
- A12_1 = +3.404e+01; A12_2 = +7.267e+01; A12_3 = +1.710e+02; A12_4 = -7.033e+02; A12_5 = +4.704e+02; A12_6 = +1.194e+01; // 2^1S -> 2^1P -- d-a
- A13_1 = +7.696e+01; A13_2 = +1.250e+02; A13_3 = +4.938e+01; A13_4 = -4.778e+02; A13_5 = +3.189e+02; A13_6 = +8.157e+00; // 2^3S -> 2^3P -- d-a
- A17_1 = +5.475e+01; A17_2 = +3.483e+05; A17_3 = +2.345e+05; A17_4 = -5.090e+05; A17_5 = +5.424e+04; A17_6 = +0.000e+00; // 2^3S -> 2^1S -- s-f
- A19_1 = +5.983e+02; A19_2 = -5.310e+02; A19_3 = +3.348e+02; A19_4 = -2.412e+02; A19_5 = +2.239e+02; A19_6 = +0.000e+00; // 2^1S -> 2^3P -- s-f
- g1 = 1; g2 = 1; g3 = 1; g4 = 1; g12 = 1; g13 = 3; g17 = 3; g19 = 1;
- printf(">> eduPIC: Setting e- / He cross sections\n");
- // load your four datafiles (make sure these names match your files!)
- CSInterpolator cs_ela ("CS/He_electron_elastic.dat"); // two‐col: E σ_ela
- CSInterpolator cs_exc1 ("CS/He_electron_exc1.dat"); // two‐col: E σ_exc1
- CSInterpolator cs_exc2 ("CS/He_electron_exc2.dat"); // two‐col: E σ_exc2
- CSInterpolator cs_ion ("CS/He_electron_ionization.dat");// two‐col: E σ_ion
- for(int i=0; i<CS_RANGES; i++){
- // your energy grid
- double en = (i==0 ? DE_CS : DE_CS * i);
- // interpolate
- sigma[E_ELA][i] = cs_ela(en);
- sigma[E_ION][i] = cs_ion(en);
- // sigma[E_EXC_1][i] = cs_exc1(en);
- // sigma[E_EXC_2][i] = cs_exc2(en);
- if (en < E_EXC_TH_1){
- sigma[E_EXC_1][i] = 0.0;
- }
- else {
- sigma[E_EXC_1][i] = forbidden_CS(A1_1, A1_2, A1_3, A1_4, A1_5, A1_6, en, E_EXC_TH_1, g1, 's');
- }
- if (en < E_EXC_TH_2){
- sigma[E_EXC_2][i] = 0.0;
- }
- else {
- sigma[E_EXC_2][i] = forbidden_CS(A2_1, A2_2, A2_3, A2_4, A2_5, A2_6, en, E_EXC_TH_2, g2, 'd');
- }
- if (en < E_EXC_TH_3){
- sigma[E_EXC_3][i] = 0.0;
- }
- else {
- sigma[E_EXC_3][i] = forbidden_CS(A3_1, A3_2, A3_3, A3_4, A3_5, A3_6, en, E_EXC_TH_3, g3, 'a');
- }
- if (en < E_EXC_TH_4){
- sigma[E_EXC_4][i] = 0.0;
- }
- else {
- sigma[E_EXC_4][i] = forbidden_CS(A4_1, A4_2, A4_3, A4_4, A4_5, A4_6, en, E_EXC_TH_4, g4, 's');
- }
- if (en < E_EXC_TH_12){
- sigma[E_EXC_12][i] = 0.0;
- }
- else {
- sigma[E_EXC_12][i] = forbidden_CS(A12_1, A12_2, A12_3, A12_4, A12_5, A12_6, en, E_EXC_TH_12, g12, 'a');
- }
- if (en < E_EXC_TH_13){
- sigma[E_EXC_13][i] = 0.0;
- }
- else {
- sigma[E_EXC_13][i] = forbidden_CS(A13_1, A13_2, A13_3, A13_4, A13_5, A13_6, en, E_EXC_TH_13, g13, 'a');
- }
- if (en < E_EXC_TH_17){
- sigma[E_EXC_17][i] = 0.0;
- }
- else {
- sigma[E_EXC_17][i] = forbidden_CS(A17_1, A17_2, A17_3, A17_4, A17_5, A17_6, en, E_EXC_TH_17, g17, 's');
- }
- if (en < E_EXC_TH_19){
- sigma[E_EXC_19][i] = 0.0;
- }
- else {
- sigma[E_EXC_19][i] = forbidden_CS(A19_1, A19_2, A19_3, A19_4, A19_5, A19_6, en, E_EXC_TH_19, g19, 's');
- }
- }
- for(int i=0; i<CS_RANGES; i++){
- // your energy grid
- double en = (i==0 ? DE_CS : DE_CS * i);
- // Superelastic for 1^1S -> 2^3S (E_SUPER_1)
- double en_plus_1 = en + E_EXC_TH_1;
- int idx1 = en_plus_1 / DE_CS;
- if (idx1 < CS_RANGES && en > 0) {
- sigma[E_SUPER_1][i] = (1.0 / 3.0) * (en_plus_1 / en) * sigma[E_EXC_1][idx1];
- } else {
- sigma[E_SUPER_1][i] = 0.0;
- }
- // Superelastic for 1^1S -> 2^1S (E_SUPER_2)
- double en_plus_2 = en + E_EXC_TH_2;
- int idx2 = en_plus_2 / DE_CS;
- if (idx2 < CS_RANGES && en > 0) {
- sigma[E_SUPER_2][i] = (1.0 / 1.0) * (en_plus_2 / en) * sigma[E_EXC_2][idx2];
- } else {
- sigma[E_SUPER_2][i] = 0.0;
- }
- // Superelastic for 1^1S -> 2^1P (E_SUPER_3)
- double en_plus_3 = en + E_EXC_TH_3;
- int idx3 = en_plus_3 / DE_CS;
- if (idx3 < CS_RANGES && en > 0) {
- sigma[E_SUPER_3][i] = (1.0 / 3.0) * (en_plus_3 / en) * sigma[E_EXC_3][idx3];
- } else {
- sigma[E_SUPER_3][i] = 0.0;
- }
- // Superelastic for 1^1S -> 2^3P (E_SUPER_4)
- double en_plus_4 = en + E_EXC_TH_4;
- int idx4 = en_plus_4 / DE_CS;
- if (idx4 < CS_RANGES && en > 0) {
- sigma[E_SUPER_4][i] = (1.0 / 9.0) * (en_plus_4 / en) * sigma[E_EXC_4][idx4];
- } else {
- sigma[E_SUPER_4][i] = 0.0;
- }
- // Superelastic for 2^1S -> 2^1P (E_SUPER_12)
- double en_plus_12 = en + E_EXC_TH_12;
- int idx12 = en_plus_12 / DE_CS;
- if (idx12 < CS_RANGES && en > 0) {
- sigma[E_SUPER_12][i] = (1.0 / 3.0) * (en_plus_12 / en) * sigma[E_EXC_12][idx12];
- } else {
- sigma[E_SUPER_12][i] = 0.0;
- }
- // Superelastic for 2^3S -> 2^3P (E_SUPER_13)
- double en_plus_13 = en + E_EXC_TH_13;
- int idx13 = en_plus_13 / DE_CS;
- if (idx13 < CS_RANGES && en > 0) {
- sigma[E_SUPER_13][i] = (3.0 / 9.0) * (en_plus_13 / en) * sigma[E_EXC_13][idx13];
- } else {
- sigma[E_SUPER_13][i] = 0.0;
- }
- // Superelastic for 2^3S -> 2^1S (E_SUPER_17)
- double en_plus_17 = en + E_EXC_TH_17;
- int idx17 = en_plus_17 / DE_CS;
- if (idx17 < CS_RANGES && en > 0) {
- sigma[E_SUPER_17][i] = (3.0 / 1.0) * (en_plus_17 / en) * sigma[E_EXC_17][idx17];
- } else {
- sigma[E_SUPER_17][i] = 0.0;
- }
- // Superelastic for 2^1S -> 2^3P (E_SUPER_19)
- double en_plus_19 = en + E_EXC_TH_19;
- int idx19 = en_plus_19 / DE_CS;
- if (idx19 < CS_RANGES && en > 0) {
- sigma[E_SUPER_19][i] = (1.0 / 9.0) * (en_plus_19 / en) * sigma[E_EXC_19][idx19];
- } else {
- sigma[E_SUPER_19][i] = 0.0;
- }
- }
- }
- //------------------------------------------------------------------------------//
- // ion cross sections: A. V. Phelps, J. Appl. Phys. 76, 747 (1994) //
- //------------------------------------------------------------------------------//
- void set_ion_cross_sections_ar(void){
- int i;
- double e_com,e_lab,qmom,qback,qiso;
- printf(">> eduPIC: Setting He+ / He cross sections\n");
- for(i=0; i<CS_RANGES; i++){
- if (i == 0) {e_com = DE_CS;} else {e_com = DE_CS * i;} // ion energy in the center of mass frame of reference
- e_lab = 2.0 * e_com; // ion energy in the laboratory frame of reference
- qiso = 7.63 *pow(10,-20) * pow(e_com, -0.5);
- qback = 1.0 * pow(10,-19) * pow( (e_com/1000.0), -0.15 ) * pow( (1.0 + e_com/1000.0), -0.25 ) * pow( (1.0 + 5.0/e_com), -0.15 );
- sigma[I_ISO][i] = qiso; // cross section for He+ / He isotropic part of elastic scattering
- sigma[I_BACK][i] = qback; // cross section for He+ / He backward elastic scattering
- }
- }
- //----------------------------------------------------------------------//
- // calculation of total cross sections for electrons and ions //
- //----------------------------------------------------------------------//
- void calc_total_cross_sections(void){
- int i;
- for(i=0; i<CS_RANGES; i++){
- sigma_tot_e[i] = (sigma[E_ELA][i] + sigma[E_EXC_1][i] + sigma[E_EXC_2][i] + sigma[E_ION][i]) * GAS_DENSITY; // total macroscopic cross section of electrons
- sigma_tot_i[i] = (sigma[I_ISO][i] + sigma[I_BACK][i]) * GAS_DENSITY; // total macroscopic cross section of ions
- }
- }
- //----------------------------------------------------------------------//
- // test of cross sections for electrons and ions //
- //----------------------------------------------------------------------//
- void test_cross_sections(void){
- FILE * f;
- int i,j;
- f = fopen("cross_sections.dat","w"); // cross sections saved in data file: cross_sections.dat
- for(i=0; i<CS_RANGES; i++){
- fprintf(f,"%12.4f ",i*DE_CS);
- for(j=0; j<N_CS; j++) fprintf(f,"%14e ",sigma[j][i]);
- fprintf(f,"\n");
- }
- fclose(f);
- }
- //---------------------------------------------------------------------//
- // find upper limit of collision frequencies //
- //---------------------------------------------------------------------//
- double max_electron_coll_freq (void){
- int i;
- double e,v,nu,nu_max;
- nu_max = 0;
- for(i=0; i<CS_RANGES; i++){
- e = i * DE_CS;
- v = sqrt(2.0 * e * EV_TO_J / E_MASS);
- nu = v * sigma_tot_e[i];
- if (nu > nu_max) {nu_max = nu;}
- }
- return nu_max;
- }
- double max_ion_coll_freq (void){
- int i;
- double e,g,nu,nu_max;
- nu_max = 0;
- for(i=0; i<CS_RANGES; i++){
- e = i * DE_CS;
- g = sqrt(2.0 * e * EV_TO_J / MU_HEHE);
- nu = g * sigma_tot_i[i];
- if (nu > nu_max) nu_max = nu;
- }
- return nu_max;
- }
- //----------------------------------------------------------------------//
- // initialization of the simulation by placing a given number of //
- // electrons and ions at random positions between the electrodes //
- //----------------------------------------------------------------------//
- // initialization of excited states distribtuion, assuming Maxwellian-Boltzmann balance -- Artem
- std::pair<double, double> init_excited_distr() {
- double part_ground = 1.0*exp(-0.0/T_neutral); // partition function for ground state
- double part_triplet = 3.0*exp(-E_EXC_TH_1*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for triplet excited state
- double part_singlet = 1.0*exp(-E_EXC_TH_2*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for singlet excited state
- double part_func_total = part_ground + part_triplet + part_singlet; // total partition function
- double n_triplet = ((part_triplet/part_func_total)*GAS_DENSITY); // denisty population of tripet state
- double n_singlet = ((part_singlet/part_func_total)*GAS_DENSITY); // density population of singlet state
- return {n_triplet, n_singlet};
- }
- void print_excitation_densities(void) {
- double total_e1 = 0.0, total_e2 = 0.0, total_e3 = 0.0, total_e4 = 0.0;
- const double cell_volume = ELECTRODE_AREA * DX;
- // Sum densities across all grid cells
- for (int p = 0; p < N_G; p++) {
- total_e1 += e1_density[p]; // 2^3S state density
- total_e2 += e2_density[p]; // 2^1S state density
- total_e3 += e3_density[p]; // 2^1P state density
- total_e4 += e4_density[p]; // 2^3P state density
- }
- std::cout << "average densities of excited states: " <<"\n";
- printf("2^3S = %12.6e | 2^1S = %12.6e\n", total_e1/N_G, total_e2/N_G);
- printf("2^1P = %12.6e | 2^3P = %12.6e\n", total_e3/N_G, total_e4/N_G);
- }
- void init(int nseed){
- int i;
- for (i=0; i<nseed; i++){
- x_e[i] = L * R01(MTgen); // initial random position of the electron
- vx_e[i] = RMB_e(MTgen); vy_e[i] = RMB_e(MTgen); vz_e[i] = RMB_e(MTgen); // initial velocity components of the electron
- x_i[i] = L * R01(MTgen); // initial random position of the ion
- vx_i[i] = RMB_n(MTgen); vy_i[i] = RMB_n(MTgen); vz_i[i] = RMB_n(MTgen); // initial velocity components of the ion
- }
- N_e = nseed; // initial number of electrons
- N_i = nseed; // initial number of ions
- // auto exc = init_excited_distr();
- // for (int p = 0; p < N_G; p++) {
- // e1_density[p] = exc.first;
- // e2_density[p] = exc.second;
- // }
- double n0 = 2.7e+9; double a0 = -1.3e+8; double b0 = 5.0e+15;
- for (int p = 0; p < N_G; p++){
- // e1_density[p] = abs(n0 * (a0 *(pow(L,4) - pow((p*DX - 0.5*L),2)) +
- // b0 * (pow(L,8) - pow((p*DX - 0.5*L),8)) ));
- // e2_density[p] = 0.5 * abs( n0 * (a0 *(pow(L,4) - pow((p*DX - 0.5*L),2)) +
- // b0 * (pow(L,8) - pow((p*DX - 0.5*L),8)) ));
- e1_density[p] = 1.0e+15;
- e2_density[p] = 1.0e+15;
- e3_density[p] = 1.0e+12; // initial non-zero safe density for 2^1P
- e4_density[p] = 1.0e+12; // initial non-zero safe density for 2^3P
- }
- }
- //----------------------------------------------------------------------//
- // e / He collision (cold gas approximation) //
- //----------------------------------------------------------------------//
- void collision_electron (double xe, double *vxe, double *vye, double *vze, int eindex){
- const double F1 = E_MASS / (E_MASS + HE_MASS);
- const double F2 = HE_MASS / (E_MASS + HE_MASS);
- double t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,rnd;
- double g,g2,gx,gy,gz,wx,wy,wz,theta,phi;
- double chi,eta,chi2,eta2,sc,cc,se,ce,st,ct,sp,cp,energy,e_sc,e_ej;
- // - Artem
- // Determine cell p where the electron is
- double c0 = xe * INV_DX;
- int p = std::max(0, std::min(N_G - 1, static_cast<int>(c0)));
- // Local densities -- Artem
- double e1_dens = e1_density[p]; // 2^3S density
- double e2_dens = e2_density[p]; // 2^1S density
- double e3_dens = e3_density[p]; // 2^1P density
- double e4_dens = e4_density[p]; // 2^3P density
- double ground_dens = GAS_DENSITY - e1_dens - e2_dens - e3_dens - e4_dens; // should I extract it? let's see Artem
- // calculate relative velocity before collision & velocity of the centre of mass
- gx = (*vxe);
- gy = (*vye);
- gz = (*vze);
- g = sqrt(gx * gx + gy * gy + gz * gz);
- wx = F1 * (*vxe);
- wy = F1 * (*vye);
- wz = F1 * (*vze);
- // find Euler angles
- if (gx == 0) {theta = 0.5 * PI;}
- else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}
- if (gy == 0) {
- if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}
- } else {phi = atan2(gz, gy);}
- st = sin(theta);
- ct = cos(theta);
- sp = sin(phi);
- cp = cos(phi);
- // choose the type of collision based on the cross sections
- // take into account energy loss in inelastic collisions
- // generate scattering and azimuth angles
- // in case of ionization handle the 'new' electron
- t0 = sigma[E_ELA ][eindex] * ground_dens;
- t1 = t0 + sigma[E_EXC_1 ][eindex] * ground_dens; // 1^1S -> 2^3S
- t2 = t1 + sigma[E_EXC_2 ][eindex] * ground_dens; // 1^1S -> 2^1S
- t3 = t2 + sigma[E_ION ][eindex] * ground_dens;
- t4 = t3 + sigma[E_SUPER_1 ][eindex] * e1_dens; // 1^1S -> 2^3S superelastic
- t5 = t4 + sigma[E_SUPER_2 ][eindex] * e2_dens; // 1^1S -> 2^1S superelastic
- t6 = t5 + sigma[E_EXC_3 ][eindex] * ground_dens; // 1^1S -> 2^1P excitation
- t7 = t6 + sigma[E_SUPER_3 ][eindex] * e3_dens; // 1^1S -> 2^1P superelastic
- t8 = t7 + sigma[E_EXC_4 ][eindex] * ground_dens; // 1^1S -> 2^3P excitation
- t9 = t8 + sigma[E_SUPER_4 ][eindex] * e4_dens; // 1^1S -> 2^3P superelastic
- t10 = t9 + sigma[E_EXC_12 ][eindex] * e2_dens; // 2^1S -> 2^1P excitation
- t11 = t10 + sigma[E_SUPER_12][eindex] * e3_dens; // 2^1S -> 2^1P superelastic
- t12 = t11 + sigma[E_EXC_13 ][eindex] * e1_dens; // 2^3S -> 2^3P excitation
- t13 = t12 + sigma[E_SUPER_13][eindex] * e4_dens; // 2^3S -> 2^3P superelastic
- t14 = t13 + sigma[E_EXC_17 ][eindex] * e1_dens; // 2^3S -> 2^1S excitation
- t15 = t14 + sigma[E_SUPER_17][eindex] * e2_dens; // 2^3S -> 2^1S superelastic
- t16 = t15 + sigma[E_EXC_19 ][eindex] * e2_dens; // 2^1S -> 2^3P excitation
- t17 = t16 + sigma[E_SUPER_19][eindex] * e4_dens; // 2^1S -> 2^3P superelastic
- rnd = R01(MTgen);
- if (rnd < (t0/t17)){ // elastic scattering
- chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering
- eta = TWO_PI * R01(MTgen); // azimuthal angle
- } else if (rnd < (t1/t17)){ // excitation 1 (triplet)
- energy = 0.5 * E_MASS * g * g; // electron energy
- energy = fabs(energy - E_EXC_TH_1 * EV_TO_J); // subtract energy loss for excitation
- g = sqrt(2.0 * energy / E_MASS); // relative velocity after energy loss
- chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering
- eta = TWO_PI * R01(MTgen); // azimuthal angle
- } else if (rnd < (t2/t17)){ // excitation 2 (singlet)
- energy = 0.5 * E_MASS * g * g; // electron energy
- energy = fabs(energy - E_EXC_TH_2 * EV_TO_J); // subtract energy loss for excitation
- g = sqrt(2.0 * energy / E_MASS); // relative velocity after energy loss
- chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering
- eta = TWO_PI * R01(MTgen); // azimuthal angle
- } else if (rnd < (t3/t17)) { // ionization
- energy = 0.5 * E_MASS * g * g; // electron energy
- energy = fabs(energy - E_ION_TH * EV_TO_J); // subtract energy loss of ionization
- // e_ej = 10.0 * tan(R01(MTgen) * atan(energy/EV_TO_J / 20.0)) * EV_TO_J; // energy of the ejected electron
- // e_sc = fabs(energy - e_ej); // energy of scattered electron after the collision
- e_ej = 0.5*energy; // energy of the ejected electron
- e_sc = fabs(energy - e_ej); // energy of scattered electron after the collision
- g = sqrt(2.0 * e_sc / E_MASS); // relative velocity of scattered electron
- g2 = sqrt(2.0 * e_ej / E_MASS); // relative velocity of ejected electron
- chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering for scattered electron (as in Turner's case)
- chi2 = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering for ejected electrons (as in Turner's case)
- eta = TWO_PI * R01(MTgen); // azimuthal angle for scattered electron
- eta2 = eta + PI; // azimuthal angle for ejected electron
- sc = sin(chi2);
- cc = cos(chi2);
- se = sin(eta2);
- ce = cos(eta2);
- gx = g2 * (ct * cc - st * sc * ce);
- gy = g2 * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
- gz = g2 * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
- x_e[N_e] = xe; // add new electron
- vx_e[N_e] = wx + F2 * gx;
- vy_e[N_e] = wy + F2 * gy;
- vz_e[N_e] = wz + F2 * gz;
- N_e++;
- x_i[N_i] = xe; // add new ion
- vx_i[N_i] = RMB_n(MTgen); // velocity is sampled from background thermal distribution
- vy_i[N_i] = RMB_n(MTgen);
- vz_i[N_i] = RMB_n(MTgen);
- N_i++;
- }
- else if (rnd < (t4/t17)) { // accounting for superelastic collisions - 1^1S -> 2^3S - Artem
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_1 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t5/t17)){ // account for superelastic collision - 1^1S -> 2^1S - Artem
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_2 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t6/t17)) { // excitation 1^1S -> 2^1P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy - E_EXC_TH_3 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t7/t17)) { // superelastic 1^1S -> 2^1P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_3 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t8/t17)) { // excitation 1^1S -> 2^3P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy - E_EXC_TH_4 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t9/t17)) { // superelastic 1^1S -> 2^3P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_4 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t10/t17)) { // excitation 2^1S -> 2^1P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy - E_EXC_TH_12 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t11/t17)) { // superelastic 2^1S -> 2^1P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_12 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t12/t17)) { // excitation 2^3S -> 2^3P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy - E_EXC_TH_13 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t13/t17)) { // superelastic 2^3S -> 2^3P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_13 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t14/t17)) { // excitation 2^3S -> 2^1S
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy - E_EXC_TH_17 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t15/t17)) { // superelastic 2^3S -> 2^1S
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_17 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else if (rnd < (t16/t17)) { // excitation 2^1S -> 2^3P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy - E_EXC_TH_19 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- } else { // superelastic 2^1S -> 2^3P
- energy = 0.5 * E_MASS * g * g;
- energy = fabs(energy + E_EXC_TH_19 * EV_TO_J);
- g = sqrt(2.0 * energy / E_MASS);
- chi = acos(1.0 - 2.0 * R01(MTgen));
- eta = TWO_PI * R01(MTgen);
- }
- // scatter the primary electron
- sc = sin(chi);
- cc = cos(chi);
- se = sin(eta);
- ce = cos(eta);
- // compute new relative velocity:
- gx = g * (ct * cc - st * sc * ce);
- gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
- gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
- // post-collision velocity of the colliding electron
- (*vxe) = wx + F2 * gx;
- (*vye) = wy + F2 * gy;
- (*vze) = wz + F2 * gz;
- }
- //----------------------------------------------------------------------//
- // He+ / He collision //
- //----------------------------------------------------------------------//
- void collision_ion (double *vx_1, double *vy_1, double *vz_1,
- double *vx_2, double *vy_2, double *vz_2, int e_index){
- double g,gx,gy,gz,wx,wy,wz,rnd;
- double theta,phi,chi,eta,st,ct,sp,cp,sc,cc,se,ce,t1,t2;
- // calculate relative velocity before collision
- // random Maxwellian target atom already selected (vx_2,vy_2,vz_2 velocity components of target atom come with the call)
- gx = (*vx_1)-(*vx_2);
- gy = (*vy_1)-(*vy_2);
- gz = (*vz_1)-(*vz_2);
- g = sqrt(gx * gx + gy * gy + gz * gz);
- wx = 0.5 * ((*vx_1) + (*vx_2));
- wy = 0.5 * ((*vy_1) + (*vy_2));
- wz = 0.5 * ((*vz_1) + (*vz_2));
- // find Euler angles
- if (gx == 0) {theta = 0.5 * PI;} else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}
- if (gy == 0) {
- if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}
- } else {phi = atan2(gz, gy);}
- // determine the type of collision based on cross sections and generate scattering angle
- t1 = sigma[I_ISO][e_index];
- t2 = t1 + sigma[I_BACK][e_index];
- rnd = R01(MTgen);
- if (rnd < (t1 /t2)){ // isotropic scattering
- chi = acos(1.0 - 2.0 * R01(MTgen)); // scattering angle
- } else { // backward scattering
- chi = PI; // scattering angle
- }
- eta = TWO_PI * R01(MTgen); // azimuthal angle
- sc = sin(chi);
- cc = cos(chi);
- se = sin(eta);
- ce = cos(eta);
- st = sin(theta);
- ct = cos(theta);
- sp = sin(phi);
- cp = cos(phi);
- // compute new relative velocity
- gx = g * (ct * cc - st * sc * ce);
- gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
- gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
- // post-collision velocity of the ion
- (*vx_1) = wx + 0.5 * gx;
- (*vy_1) = wy + 0.5 * gy;
- (*vz_1) = wz + 0.5 * gz;
- }
- //-----------------------------------------------------------------//
- // solve Poisson equation (Thomas algorithm) //
- //-----------------------------------------------------------------//
- void solve_Poisson (xvector rho1, double tt){
- const double A = 1.0;
- const double B = -2.0;
- const double C = 1.0;
- const double S = 1.0 / (2.0 * DX);
- const double ALPHA = -DX * DX / EPSILON0;
- xvector g, w, f;
- int i;
- // apply potential to the electrodes - boundary conditions
- pot[0] = VOLTAGE * cos(OMEGA * tt); // potential at the powered electrode
- pot[N_G-1] = 0.0; // potential at the grounded electrode
- // solve Poisson equation
- for(i=1; i<=N_G-2; i++) f[i] = ALPHA * rho1[i];
- f[1] -= pot[0];
- f[N_G-2] -= pot[N_G-1];
- w[1] = C/B;
- g[1] = f[1]/B;
- for(i=2; i<=N_G-2; i++){
- w[i] = C / (B - A * w[i-1]);
- g[i] = (f[i] - A * g[i-1]) / (B - A * w[i-1]);
- }
- pot[N_G-2] = g[N_G-2];
- for (i=N_G-3; i>0; i--) pot[i] = g[i] - w[i] * pot[i+1]; // potential at the grid points between the electrodes
- // compute electric field
- for(i=1; i<=N_G-2; i++) efield[i] = (pot[i-1] - pot[i+1]) * S; // electric field at the grid points between the electrodes
- efield[0] = (pot[0] - pot[1]) * INV_DX - rho1[0] * DX / (2.0 * EPSILON0); // powered electrode
- efield[N_G-1] = (pot[N_G-2] - pot[N_G-1]) * INV_DX + rho1[N_G-1] * DX / (2.0 * EPSILON0); // grounded electrode
- }
- //---------------------------------------------------------------------//
- // simulation of one radiofrequency cycle //
- //---------------------------------------------------------------------//
- void accumulate_rates() {
- double v_sqr, velocity, energy, c0_temp;
- int energy_index, p_temp;
- const double Volume = (ELECTRODE_AREA * DX);
- for (int k=0; k<N_e; k++){
- v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
- velocity = sqrt(v_sqr);
- energy = 0.5 * E_MASS * v_sqr / EV_TO_J;
- energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
- c0_temp = x_e[k] * INV_DX;
- p_temp = std::max(0, std::min(N_G-1, static_cast<int>(c0_temp)));
- sum_electron_density[p_temp] += WEIGHT / Volume;
- // 1^1S -> 2^3S
- sum_rate1f[p_temp] += sigma[E_EXC_1][energy_index] * velocity * WEIGHT;
- sum_rate1b[p_temp] += sigma[E_SUPER_1][energy_index] * velocity * WEIGHT;
- // 1^1S -> 2^1S
- sum_rate2f[p_temp] += sigma[E_EXC_2][energy_index] * velocity * WEIGHT;
- sum_rate2b[p_temp] += sigma[E_SUPER_2][energy_index] * velocity * WEIGHT;
- // 1^1S -> 2^1P
- sum_rate3f[p_temp] += sigma[E_EXC_3][energy_index] * velocity * WEIGHT;
- sum_rate3b[p_temp] += sigma[E_SUPER_3][energy_index] * velocity * WEIGHT;
- // 1^1S -> 2^3P
- sum_rate4f[p_temp] += sigma[E_EXC_4][energy_index] * velocity * WEIGHT;
- sum_rate4b[p_temp] += sigma[E_SUPER_4][energy_index] * velocity * WEIGHT;
- // 2^1S -> 2^1P
- sum_rate12f[p_temp] += sigma[E_EXC_12][energy_index] * velocity * WEIGHT;
- sum_rate12b[p_temp] += sigma[E_SUPER_12][energy_index] * velocity * WEIGHT;
- // 2^3S -> 2^3P
- sum_rate13f[p_temp] += sigma[E_EXC_13][energy_index] * velocity * WEIGHT;
- sum_rate13b[p_temp] += sigma[E_SUPER_13][energy_index] * velocity * WEIGHT;
- // 2^3S -> 2^1S
- sum_rate17f[p_temp] += sigma[E_EXC_17][energy_index] * velocity * WEIGHT;
- sum_rate17b[p_temp] += sigma[E_SUPER_17][energy_index] * velocity * WEIGHT;
- // 2^1S -> 2^3P
- sum_rate19f[p_temp] += sigma[E_EXC_19][energy_index] * velocity * WEIGHT;
- sum_rate19b[p_temp] += sigma[E_SUPER_19][energy_index] * velocity * WEIGHT;
- }
- }
- // averaging the rates each RF cycle
- void average_rates() {
- const double inv_NT = 1.0 / (N_avg * N_T); // averaging through N_avg RF periods each contains N_T cycles
- const double inv_Volume = 1.0/(ELECTRODE_AREA * DX);
- for (int p = 0; p < N_G; p++){
- avg_electron_density[p] = sum_electron_density[p] * inv_NT;
- }
- for (int p = 0; p < N_G; p++) {
- avg_rate1f[p] = sum_rate1f[p] * inv_NT * inv_Volume;
- avg_rate1b[p] = sum_rate1b[p] * inv_NT * inv_Volume;
- avg_rate2f[p] = sum_rate2f[p] * inv_NT * inv_Volume;
- avg_rate2b[p] = sum_rate2b[p] * inv_NT * inv_Volume;
- avg_rate3f[p] = sum_rate3f[p] * inv_NT * inv_Volume;
- avg_rate3b[p] = sum_rate3b[p] * inv_NT * inv_Volume;
- avg_rate4f[p] = sum_rate4f[p] * inv_NT * inv_Volume;
- avg_rate4b[p] = sum_rate4b[p] * inv_NT * inv_Volume;
- avg_rate12f[p] = sum_rate12f[p] * inv_NT * inv_Volume;
- avg_rate12b[p] = sum_rate12b[p] * inv_NT * inv_Volume;
- avg_rate13f[p] = sum_rate13f[p] * inv_NT * inv_Volume;
- avg_rate13b[p] = sum_rate13b[p] * inv_NT * inv_Volume;
- avg_rate17f[p] = sum_rate17f[p] * inv_NT * inv_Volume;
- avg_rate17b[p] = sum_rate17b[p] * inv_NT * inv_Volume;
- avg_rate19f[p] = sum_rate19f[p] * inv_NT * inv_Volume;
- avg_rate19b[p] = sum_rate19b[p] * inv_NT * inv_Volume;
- }
- }
- void update_excited_dens() {
- xvector e1_density_bw, e2_density_bw, e3_density_bw, e4_density_bw, gas_dens_local;
- double S0, S1, S2, S3, S4; // source terms to account for reaction
- double diff1, diff2; // diffusion terms only for metastable states
- const double v_mean = sqrt(8.0*K_BOLTZMANN*T_neutral/(M_PI*HE_MASS)); // mean spead of He excited states, assuming T~Tg
- double k20f, k21f, k26f;
- k20f = 1.799e+09; k21f = 1.976e+06; k26f = 1.022e+07; // radiative transitions
- // k20f = 0.000e+00; k21f = 0.000e+00; k26f = 0.000e+00; // radiative transitions
- bool converged = false;
- double accuracy = 1.0E-6; // demanded accuracy
- double alpha = 0.01;
- double dt = alpha*DX*DX/Diff_coeff;
- int max_iter = 1E8;
- // constants to calculate reabsorption:
- double C1 = sqrt(2.0*K_BOLTZMANN*T_neutral/HE_MASS);
- double C2 = E_CHARGE*E_CHARGE/(4.0*EPSILON0*E_MASS*C_light);
- // spectral line 2^1P -> 1^1 S
- double lambda20, d_nu20, phi20, f_lu20, alpha20, tau20, absorb20; // spectral line constants for for spectral line 2^1P -> 1^1S
- lambda20 = 584.334357e-10; // wavelength of transition
- d_nu20 = 1.0/lambda20 * C1; // Doppler broadening delta_nu
- phi20 = 1.0/(d_nu20*sqrt(M_PI)); // Doppler broadening phi
- f_lu20 = 2.7616e-01; // oscillator strength
- double lambda21, d_nu21, phi21, f_lu21, alpha21, tau21, absorb21; // spectral line constants for for spectral line 2^1P -> 2^1S
- lambda21 = 20587.0e-10;
- d_nu21 = 1.0/lambda21 * C1;
- phi21 = 1.0/(d_nu21*sqrt(M_PI));
- f_lu21 = 3.764e-01;
- double lambda26, d_nu26, phi26, f_lu26, alpha26, tau26, absorb26; // spectral line constants for for spectral line 2^3P -> 2^3S
- lambda26 = 10830.0e-10;
- d_nu26 = 1.0/lambda26 * C1;
- phi26 = 1.0/(d_nu26*sqrt(M_PI));
- f_lu26 = 5.3907e-01;
- // creating a helping array for backward step and also works with initial condition
- std::memcpy(e1_density_bw, e1_density, sizeof(e1_density));
- std::memcpy(e2_density_bw, e2_density, sizeof(e2_density));
- std::memcpy(e3_density_bw, e3_density, sizeof(e3_density));
- std::memcpy(e4_density_bw, e4_density, sizeof(e4_density));
- for (int i = 0; i < N_G; i++){
- gas_dens_local[i] = GAS_DENSITY;
- }
- double rcmb_coeff = gamma_k*v_mean/(2.0*(2.0-gamma_k));
- for (int idx = 0; idx < max_iter; idx++) {
- // swapping the arrays
- std::swap(e1_density_bw, e1_density);
- std::swap(e2_density_bw, e2_density);
- std::swap(e3_density_bw, e3_density);
- std::swap(e4_density_bw, e4_density);
- for (int i = 1; i < N_G-1; i++){
- // calculate source terms of diffusion equation, S0 for neutral atoms of He
- // S1 for triplet excited state, S2 for singlet excited state
- // triplet and singlet are diffusing, neutral atoms are not
- // S0 = dt*(-(avg_rate1f[i]+avg_rate2f[i] + avg_rate3f[i] + avg_rate4f[i])*gas_dens_local[i]
- // + avg_rate1b[i]*e1_density_bw[i] + avg_rate2b[i]*e2_density_bw[i]
- // + avg_rate3b[i]*e3_density_bw[i] + avg_rate4b[i]*e4_density_bw[i] + k20f*e3_density_bw[i]);
- // calculating reabsorption coefficient for 2^1P -> 1^1S
- alpha20 = C2 * f_lu20 * gas_dens_local[i] * (1.0 - (e3_density_bw[i]*1)/(gas_dens_local[i]*3)) * phi20;
- tau20 = alpha20*L/2.0;
- absorb20 = compute_escape_factor(tau20);
- // calculating reabsorption coefficient for 2^1P -> 2^1S
- alpha21 = C2 * f_lu21 * e2_density_bw[i] * (1.0 - (e3_density_bw[i]*1)/(e2_density_bw[i]*3)) * phi21;
- tau21 = alpha21*L/2.0;
- absorb21 = compute_escape_factor(tau21);
- // calculating reabsorption coefficient for 2^3P -> 2^3S
- alpha26 = C2 * f_lu26 * e1_density_bw[i] * (1.0 - (e4_density_bw[i]*3)/(e1_density_bw[i]*9)) * phi26;
- tau26 = alpha26*L/2.0;
- absorb26 = compute_escape_factor(tau26);
- S1 = dt*( avg_rate1f[i]*gas_dens_local[i] + avg_rate13b[i]*e4_density_bw[i] + avg_rate17b[i]*e2_density_bw[i]
- - (avg_rate1b[i] + avg_rate13f[i] + avg_rate17f[i])*e1_density_bw[i] + k26f*absorb26*e4_density_bw[i] ) ;
- S2 = dt*(avg_rate2f[i]*gas_dens_local[i] + avg_rate12b[i]*e3_density_bw[i]
- + avg_rate19b[i]*e4_density_bw[i] + avg_rate17f[i]*e1_density_bw[i]
- - (avg_rate2b[i] + avg_rate12f[i] + avg_rate19f[i] + avg_rate17b[i])*e2_density_bw[i] + k21f*absorb21*e3_density_bw[i] );
- S3 = dt*(avg_rate3f[i]*gas_dens_local[i] + avg_rate12f[i]*e2_density_bw[i]
- - (avg_rate3b[i] + avg_rate12b[i])*e3_density_bw[i] - (k20f*absorb20+k21f*absorb21)*e3_density_bw[i] );
- S4 = dt*(avg_rate4f[i]*gas_dens_local[i] + avg_rate13f[i]*e1_density_bw[i] + avg_rate19f[i]*e2_density_bw[i]
- - (avg_rate4b[i] + avg_rate13b[i] + avg_rate19b[i])*e4_density_bw[i] - k26f*absorb26*e4_density_bw[i]);
- diff1 = e1_density_bw[i]*(1.0 - 2.0*alpha) + alpha*(e1_density_bw[i-1] + e1_density_bw[i+1]);
- diff2 = e2_density_bw[i]*(1.0 - 2.0*alpha) + alpha*(e2_density_bw[i-1] + e2_density_bw[i+1]);
- e1_density[i] = diff1 + S1;
- e2_density[i] = diff2 + S2;
- e3_density[i] = e3_density_bw[i] + S3;
- e4_density[i] = e4_density_bw[i] + S4;
- // gas_dens_local[i] += S0;
- // Prevent negative densities
- e1_density[i] = std::max(e1_density[i], 0.0);
- e2_density[i] = std::max(e2_density[i], 0.0);
- e3_density[i] = std::max(e3_density[i], 0.0);
- e4_density[i] = std::max(e4_density[i], 0.0);
- }
- // boundary conditions: now it's zero so no excited states density at the electrodes, yield is zero
- e1_density[0] = e1_density[1]/(1.0 + rcmb_coeff*DX/Diff_coeff);
- e1_density[N_G-1] = e1_density[N_G-2]/(1.0 + rcmb_coeff*DX/Diff_coeff);
- e2_density[0] = e2_density[1]/(1.0 + rcmb_coeff*DX/Diff_coeff);
- e2_density[N_G-1] = e2_density[N_G-2]/(1.0 + rcmb_coeff*DX/Diff_coeff);
- e3_density[0] = e3_density_bw[0] + dt*(avg_rate3f[0]*gas_dens_local[0] + avg_rate12f[0]*e2_density_bw[0]
- - (avg_rate3b[0] + avg_rate12b[0])*e3_density_bw[0] - (k20f+k21f)*e3_density_bw[0] );
- e3_density[N_G-1] = e3_density_bw[N_G-1] + dt*(avg_rate3f[N_G-1]*gas_dens_local[N_G-1] + avg_rate12f[N_G-1]*e2_density_bw[N_G-1]
- - (avg_rate3b[N_G-1] + avg_rate12b[N_G-1])*e3_density_bw[N_G-1] - (k20f+k21f)*e3_density_bw[N_G-1] );
- e4_density[0] = e4_density_bw[0] + dt*(avg_rate4f[0]*gas_dens_local[0] + avg_rate13f[0]*e1_density_bw[0] + avg_rate19f[0]*e2_density_bw[0]
- - (avg_rate4b[0] + avg_rate13b[0] + avg_rate19b[0])*e4_density_bw[0] - k26f*e4_density_bw[0]);
- e4_density[N_G-1] = e4_density_bw[N_G-1] + dt*(avg_rate4f[N_G-1]*gas_dens_local[N_G-1] + avg_rate13f[N_G-1]*e1_density_bw[N_G-1] + avg_rate19f[N_G-1]*e2_density_bw[N_G-1]
- - (avg_rate4b[N_G-1] + avg_rate13b[N_G-1] + avg_rate19b[N_G-1])*e4_density_bw[N_G-1] - k26f*e4_density_bw[N_G-1]);
- //checking convergence
- double rel1 = 0.0, rel2 = 0.0, rel3 = 0.0, rel4 = 0.0;
- for (int i = 1; i < N_G-1; ++i) {
- double abs1 = fabs(e1_density[i] - e1_density_bw[i]);
- double abs2 = fabs(e2_density[i] - e2_density_bw[i]);
- double abs3 = fabs(e3_density[i] - e3_density_bw[i]);
- double abs4 = fabs(e4_density[i] - e4_density_bw[i]);
- double base1 = std::max(e1_density[i], 1e-12);
- double base2 = std::max(e2_density[i], 1e-12);
- double base3 = std::max(e3_density[i], 1e-12);
- double base4 = std::max(e4_density[i], 1e-12);
- rel1 = std::max(rel1, abs1 / e1_density_bw[i]);
- rel2 = std::max(rel2, abs2 / e2_density_bw[i]);
- rel3 = std::max(rel3, abs3 / e3_density_bw[i]);
- rel4 = std::max(rel4, abs4 / e4_density_bw[i]);
- }
- if (rel1 < accuracy && rel2 < accuracy && rel3 < accuracy && rel4 < accuracy) {
- converged = true;
- std::cout << "Steady-state reached after " << idx << " iterations.\n";
- std::cout << "Relative changes are: " << rel1 << " " << rel2 << " " << rel3 << " " << rel4 << "\n";
- break;
- }
- }
- if (!converged) {
- std::cerr << "Steady-state not reached after " << max_iter << " iterations.\n";
- }
- }
- void do_one_cycle (void){
- const double DV = ELECTRODE_AREA * DX;
- const double FACTOR_W = WEIGHT / DV;
- const double FACTOR_E = DT_E / E_MASS * E_CHARGE;
- const double FACTOR_I = DT_I / HE_MASS * E_CHARGE;
- const double MIN_X = 0.45 * L; // min. position for EEPF collection
- const double MAX_X = 0.55 * L; // max. position for EEPF collection
- int k, t, p, energy_index;
- double g, g_sqr, gx, gy, gz, vx_a, vy_a, vz_a, e_x, energy, nu, p_coll, v_sqr, velocity;
- double mean_v, c0, c1, c2, rate;
- bool out;
- xvector rho;
- int t_index;
- nu_avg = 0.0;
- for (t=0; t<N_T; t++){ // the RF period is divided into N_T equal time intervals (time step DT_E)
- Time += DT_E; // update of the total simulated time
- t_index = t / N_BIN; // index for XT distributions
- // step 1: compute densities at grid points
- for(p=0; p<N_G; p++) e_density[p] = 0; // electron density - computed in every time step
- for(k=0; k<N_e; k++){
- if (p < 0) p = 0;
- else if (p > N_G - 2) p = N_G - 2;
- c0 = x_e[k] * INV_DX;
- p = int(c0);
- e_density[p] += (p + 1 - c0) * FACTOR_W;
- e_density[p+1] += (c0 - p) * FACTOR_W;
- }
- e_density[0] *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)
- e_density[N_G-1] *= 2.0;
- for(p=0; p<N_G; p++) cumul_e_density[p] += e_density[p];
- if ((t % N_SUB) == 0) { // ion density - computed in every N_SUB-th time steps (subcycling)
- for(p=0; p<N_G; p++) i_density[p] = 0;
- for(k=0; k<N_i; k++){
- c0 = x_i[k] * INV_DX;
- p = int(c0);
- i_density[p] += (p + 1 - c0) * FACTOR_W;
- i_density[p+1] += (c0 - p) * FACTOR_W;
- }
- i_density[0] *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)
- i_density[N_G-1] *= 2.0;
- }
- for(p=0; p<N_G; p++) cumul_i_density[p] += i_density[p];
- // step 2: solve Poisson equation
- for(p=0; p<N_G; p++) rho[p] = E_CHARGE * (i_density[p] - e_density[p]); // get charge density
- solve_Poisson(rho,Time); // compute potential and electric field
- // steps 3 & 4: move particles according to electric field interpolated to particle positions
- for(k=0; k<N_e; k++){ // move all electrons in every time step
- c0 = x_e[k] * INV_DX;
- p = int(c0);
- c1 = p + 1.0 - c0;
- c2 = c0 - p;
- e_x = c1 * efield[p] + c2 * efield[p+1];
- if (measurement_mode) {
- // measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'
- mean_v = vx_e[k] - 0.5 * e_x * FACTOR_E;
- counter_e_xt[p][t_index] += c1;
- counter_e_xt[p+1][t_index] += c2;
- ue_xt[p][t_index] += c1 * mean_v;
- ue_xt[p+1][t_index] += c2 * mean_v;
- v_sqr = mean_v * mean_v + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
- energy = 0.5 * E_MASS * v_sqr / EV_TO_J;
- meanee_xt[p][t_index] += c1 * energy;
- meanee_xt[p+1][t_index] += c2 * energy;
- energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
- velocity = sqrt(v_sqr);
- double local_neut_dens = GAS_DENSITY - e1_density[p] - e2_density[p];
- rate = sigma[E_ION][energy_index] * velocity * DT_E * local_neut_dens;
- ioniz_rate_xt[p][t_index] += c1 * rate;
- ioniz_rate_xt[p+1][t_index] += c2 * rate;
- // measure EEPF in the center
- if ((MIN_X < x_e[k]) && (x_e[k] < MAX_X)){
- energy_index = (int)(energy / DE_EEPF);
- if (energy_index < N_EEPF) {eepf[energy_index] += 1.0;}
- mean_energy_accu_center += energy;
- mean_energy_counter_center++;
- }
- }
- // update velocity and position
- vx_e[k] -= e_x * FACTOR_E;
- x_e[k] += vx_e[k] * DT_E;
- }
- if ((t % N_SUB) == 0) { // move all ions in every N_SUB-th time steps (subcycling)
- for(k=0; k<N_i; k++){
- c0 = x_i[k] * INV_DX;
- p = int(c0);
- c1 = p + 1 - c0;
- c2 = c0 - p;
- e_x = c1 * efield[p] + c2 * efield[p+1];
- if (measurement_mode) {
- // measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'
- mean_v = vx_i[k] + 0.5 * e_x * FACTOR_I;
- counter_i_xt[p][t_index] += c1;
- counter_i_xt[p+1][t_index] += c2;
- ui_xt[p][t_index] += c1 * mean_v;
- ui_xt[p+1][t_index] += c2 * mean_v;
- v_sqr = mean_v * mean_v + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
- energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;
- meanei_xt[p][t_index] += c1 * energy;
- meanei_xt[p+1][t_index] += c2 * energy;
- }
- // update velocity and position and accumulate absorbed energy
- vx_i[k] += e_x * FACTOR_I;
- x_i[k] += vx_i[k] * DT_I;
- }
- }
- // step 5: check boundaries -- changed by Artem
- k = 0;
- while(k < N_e) { // check boundaries for all electrons in every time step
- out = false;
- if ((x_e[k] < 0) || (x_e[k] > L)) {out = true;} // the electron is out at the electrodes
- if (out) {
- if (R01(MTgen)< r) { // reflecting elastically with r probability
- x_e[k] < 0 ? x_e[k] *= -1.0 : x_e[k] = 2.0*L - x_e[k]; // if left/right border -> put on a left/right border
- vx_e[k] *= -1.0; // elastic reflection
- }
- else { // remove the electron
- x_e[k] < 0 ? N_e_abs_pow++ : N_e_abs_gnd++; // absorption calculation at the powered/grounded (left/right) electrode
- x_e [k] = x_e [N_e-1]; // pushing last element on a vacant place
- vx_e[k] = vx_e[N_e-1];
- vy_e[k] = vy_e[N_e-1];
- vz_e[k] = vz_e[N_e-1];
- N_e--;
- }
- } else k++;
- }
- if ((t % N_SUB) == 0) { // check boundaries for all ions in every N_SUB-th time steps (subcycling)
- k = 0;
- while(k < N_i) {
- out = false;
- if (x_i[k] < 0) { // the ion is out at the powered electrode
- N_i_abs_pow++;
- out = true;
- v_sqr = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
- energy = 0.5 * HE_MASS * v_sqr/ EV_TO_J;
- energy_index = (int)(energy / DE_IFED);
- if (energy_index < N_IFED) {ifed_pow[energy_index]++;} // save IFED at the powered electrode
- if (R01(MTgen)< gamma_i) {
- if (N_e < MAX_N_P){ // adding a new electron on the left electrode
- N_e++; // with energy sampled from Maxwellian distribution
- vx_e[N_e-1] = RMB_e(MTgen); vy_e[N_e-1] = RMB_e(MTgen); vz_e[N_e-1] = RMB_e(MTgen);
- x_e[N_e-1] = 0.0;
- }
- else{
- std::cerr << "Impossible to emit a new SE at the powered electrode" <<"\n";
- }
- }
- }
- if (x_i[k] > L) { // the ion is out at the grounded electrode
- N_i_abs_gnd++;
- out = true;
- v_sqr = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
- energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;
- energy_index = (int)(energy / DE_IFED);
- if (energy_index < N_IFED) {ifed_gnd[energy_index]++;} // save IFED at the grounded electrode
- if (R01(MTgen)< gamma_i) {
- if (N_e < MAX_N_P){ // adding a new electron on the left electrode
- N_e++; // with energy sampled from Maxwellian distribution
- vx_e[N_e-1] = RMB_e(MTgen); vy_e[N_e-1] = RMB_e(MTgen); vz_e[N_e-1] = RMB_e(MTgen);
- x_e[N_e-1] = L;
- }
- else{
- std::cerr << "Impossible to emit a new SE at the grounded electrode" <<"\n";
- }
- }
- }
- if (out) { // delete the ion, if out
- x_i [k] = x_i [N_i-1];
- vx_i[k] = vx_i[N_i-1];
- vy_i[k] = vy_i[N_i-1];
- vz_i[k] = vz_i[N_i-1];
- N_i--;
- } else k++;
- }
- }
- // step 6: collisions
- for (k=0; k<N_e; k++){ // checking for occurrence of a collision for all electrons in every time step
- v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
- velocity = sqrt(v_sqr);
- energy = 0.5 * E_MASS * v_sqr / EV_TO_J;
- energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
- // double ground_dens_local = GAS_DENSITY; // let the ground density be fixed in each cell according to our DRR solver
- // Artem - adding superelastic impact on total collisoinal probability//
- int e_crdnt = std::max(0, std::min(N_G-1, static_cast<int>(x_e[k] * INV_DX)));
- double sigma_super_e = sigma[E_SUPER_1][energy_index] * e1_density[e_crdnt] + sigma[E_SUPER_2][energy_index] * e2_density[e_crdnt]
- + sigma[E_SUPER_3][energy_index] * e3_density[e_crdnt] + sigma[E_SUPER_4][energy_index] * e4_density[e_crdnt]
- + sigma[E_SUPER_12][energy_index] * e3_density[e_crdnt] + sigma[E_SUPER_13][energy_index] * e4_density[e_crdnt]
- + sigma[E_SUPER_17][energy_index] * e2_density[e_crdnt] + sigma[E_SUPER_19][energy_index] * e4_density[e_crdnt];
- double ground_dens_local = GAS_DENSITY - e1_density[e_crdnt] - e2_density[e_crdnt] - e3_density[e_crdnt] - e4_density[e_crdnt];
- double sigma_ground = (sigma[E_ELA][energy_index] + sigma[E_EXC_1][energy_index]
- + sigma[E_EXC_2][energy_index] + sigma[E_ION][energy_index]
- + sigma[E_EXC_3][energy_index] + sigma[E_EXC_4][energy_index]) * ground_dens_local;
- double sigma_stepwise = sigma[E_EXC_12][energy_index] * e2_density[e_crdnt] + sigma[E_EXC_13][energy_index] * e1_density[e_crdnt]
- + sigma[E_EXC_17][energy_index] * e1_density[e_crdnt] + sigma[E_EXC_19][energy_index] * e2_density[e_crdnt];
- nu = (sigma_ground + sigma_super_e + sigma_stepwise) * velocity;
- nu_avg += nu;
- p_coll = 1 - exp(- nu * DT_E); // collision probability for electrons
- if (R01(MTgen) < p_coll) { // electron collision takes place
- collision_electron(x_e[k], &vx_e[k], &vy_e[k], &vz_e[k], energy_index);
- N_e_coll++;
- }
- }
- if ((t % N_SUB) == 0) { // checking for occurrence of a collision for all ions in every N_SUB-th time steps (subcycling)
- for (k=0; k<N_i; k++){
- vx_a = RMB_n(MTgen); // pick velocity components of a random target gas atom
- vy_a = RMB_n(MTgen);
- vz_a = RMB_n(MTgen);
- gx = vx_i[k] - vx_a; // compute the relative velocity of the collision partners
- gy = vy_i[k] - vy_a;
- gz = vz_i[k] - vz_a;
- g_sqr = gx * gx + gy * gy + gz * gz;
- g = sqrt(g_sqr);
- energy = 0.5 * MU_HEHE * g_sqr / EV_TO_J;
- energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
- nu = sigma_tot_i[energy_index] * g;
- p_coll = 1 - exp(- nu * DT_I); // collision probability for ions
- if (R01(MTgen)< p_coll) { // ion collision takes place
- collision_ion (&vx_i[k], &vy_i[k], &vz_i[k], &vx_a, &vy_a, &vz_a, energy_index);
- N_i_coll++;
- }
- }
- }
- //step 7: accumulate rates
- accumulate_rates();
- if (measurement_mode) {
- // collect 'xt' data from the grid
- for (p=0; p<N_G; p++) {
- pot_xt [p][t_index] += pot[p];
- efield_xt[p][t_index] += efield[p];
- ne_xt [p][t_index] += e_density[p];
- ni_xt [p][t_index] += i_density[p];
- e1_xt [p][t_index] += e1_density[p]; // Artem
- e2_xt [p][t_index] += e2_density[p]; // Artem
- }
- }
- if ((t % 1000) == 0) printf(" c = %8d t = %8d #e = %8d #i = %8d\n", cycle,t,N_e,N_i);
- }
- counter++;
- // updating denisites each N_avg cycles: --- Artem
- if (counter%N_avg == 0) {
- // compute average rates over a cycle
- average_rates();
- // updating densities
- update_excited_dens();
- // reset accumulators
- memset(sum_rate1f, 0, sizeof sum_rate1f);
- memset(sum_rate1b, 0, sizeof sum_rate1b);
- memset(sum_rate2f, 0, sizeof sum_rate2f);
- memset(sum_rate2b, 0, sizeof sum_rate2b);
- memset(sum_rate3f, 0, sizeof sum_rate3f);
- memset(sum_rate3b, 0, sizeof sum_rate3b);
- memset(sum_rate4f, 0, sizeof sum_rate4f);
- memset(sum_rate4b, 0, sizeof sum_rate4b);
- memset(sum_rate12f, 0, sizeof sum_rate12f);
- memset(sum_rate12b, 0, sizeof sum_rate12b);
- memset(sum_rate13f, 0, sizeof sum_rate13f);
- memset(sum_rate13b, 0, sizeof sum_rate13b);
- memset(sum_rate17f, 0, sizeof sum_rate17f);
- memset(sum_rate17b, 0, sizeof sum_rate17b);
- memset(sum_rate19f, 0, sizeof sum_rate19f);
- memset(sum_rate19b, 0, sizeof sum_rate19b);
- memset(sum_electron_density, 0, sizeof sum_electron_density);
- }
- //calculate nu electron average:
- nu_avg /= (N_T*N_e);
- fprintf(datafile,"%8d %8d %8d\n",cycle,N_e,N_i);
- print_excitation_densities();
- }
- //---------------------------------------------------------------------//
- // save particle coordinates //
- //---------------------------------------------------------------------//
- void save_particle_data(){
- double d;
- FILE * f;
- char fname[80];
- strcpy(fname,"picdata.bin");
- f = fopen(fname,"wb");
- fwrite(&Time,sizeof(double),1,f);
- d = (double)(cycles_done);
- fwrite(&d,sizeof(double),1,f);
- d = (double)(N_e);
- fwrite(&d,sizeof(double),1,f);
- fwrite(x_e, sizeof(double),N_e,f);
- fwrite(vx_e,sizeof(double),N_e,f);
- fwrite(vy_e,sizeof(double),N_e,f);
- fwrite(vz_e,sizeof(double),N_e,f);
- d = (double)(N_i);
- fwrite(&d,sizeof(double),1,f);
- fwrite(x_i, sizeof(double),N_i,f);
- fwrite(vx_i,sizeof(double),N_i,f);
- fwrite(vy_i,sizeof(double),N_i,f);
- fwrite(vz_i,sizeof(double),N_i,f);
- // saving excited state densities and rates - Artem
- fwrite(e1_density, sizeof(double), N_G, f); fwrite(e2_density, sizeof(double), N_G, f);
- fwrite(e3_density, sizeof(double), N_G, f); fwrite(e4_density, sizeof(double), N_G, f);
- fwrite(sum_rate1f, sizeof(double), N_G, f); fwrite(sum_rate1b, sizeof(double), N_G, f);
- fwrite(sum_rate2f, sizeof(double), N_G, f); fwrite(sum_rate2b, sizeof(double), N_G, f);
- fwrite(sum_rate3f, sizeof(double), N_G, f); fwrite(sum_rate3b, sizeof(double), N_G, f);
- fwrite(sum_rate4f, sizeof(double), N_G, f); fwrite(sum_rate4b, sizeof(double), N_G, f);
- fwrite(sum_rate12f, sizeof(double), N_G, f); fwrite(sum_rate12b, sizeof(double), N_G, f);
- fwrite(sum_rate13f, sizeof(double), N_G, f); fwrite(sum_rate13b, sizeof(double), N_G, f);
- fwrite(sum_rate17f, sizeof(double), N_G, f); fwrite(sum_rate17b, sizeof(double), N_G, f);
- fwrite(sum_rate19f, sizeof(double), N_G, f); fwrite(sum_rate19b, sizeof(double), N_G, f);
- fwrite(avg_rate1f, sizeof(double), N_G, f); fwrite(avg_rate1b, sizeof(double), N_G, f);
- fwrite(avg_rate2f, sizeof(double), N_G, f); fwrite(avg_rate2b, sizeof(double), N_G, f);
- fwrite(avg_rate3f, sizeof(double), N_G, f); fwrite(avg_rate3b, sizeof(double), N_G, f);
- fwrite(avg_rate4f, sizeof(double), N_G, f); fwrite(avg_rate4b, sizeof(double), N_G, f);
- fwrite(avg_rate12f, sizeof(double), N_G, f); fwrite(avg_rate12b, sizeof(double), N_G, f);
- fwrite(avg_rate13f, sizeof(double), N_G, f); fwrite(avg_rate13b, sizeof(double), N_G, f);
- fwrite(avg_rate17f, sizeof(double), N_G, f); fwrite(avg_rate17b, sizeof(double), N_G, f);
- fwrite(avg_rate19f, sizeof(double), N_G, f); fwrite(avg_rate19b, sizeof(double), N_G, f);
- fwrite(sum_electron_density, sizeof(double), N_G, f);
- fwrite(avg_electron_density, sizeof(double), N_G, f);
- fclose(f);
- printf(">> eduPIC: data saved : %d electrons %d ions, excited states densities , %d cycles completed, time is %e [s]\n",N_e,N_i,cycles_done,Time);
- }
- //---------------------------------------------------------------------//
- // load particle coordinates //
- //---------------------------------------------------------------------//
- void load_particle_data(){
- double d;
- FILE * f;
- char fname[80];
- strcpy(fname,"picdata.bin");
- f = fopen(fname,"rb");
- if (f==NULL) {printf(">> eduPIC: ERROR: No particle data file found, try running initial cycle using argument '0'\n"); exit(0); }
- fread(&Time,sizeof(double),1,f);
- fread(&d,sizeof(double),1,f);
- cycles_done = int(d);
- fread(&d,sizeof(double),1,f);
- N_e = int(d);
- fread(x_e, sizeof(double),N_e,f);
- fread(vx_e,sizeof(double),N_e,f);
- fread(vy_e,sizeof(double),N_e,f);
- fread(vz_e,sizeof(double),N_e,f);
- fread(&d,sizeof(double),1,f);
- N_i = int(d);
- fread(x_i, sizeof(double),N_i,f);
- fread(vx_i,sizeof(double),N_i,f);
- fread(vy_i,sizeof(double),N_i,f);
- fread(vz_i,sizeof(double),N_i,f);
- // reading excited states densities and rates info -- Artem
- fread(e1_density, sizeof(double), N_G, f);
- fread(e2_density, sizeof(double), N_G, f);
- fread(e3_density, sizeof(double), N_G, f);
- fread(e4_density, sizeof(double), N_G, f);
- fread(sum_rate1f, sizeof(double), N_G, f); fread(sum_rate1b, sizeof(double), N_G, f);
- fread(sum_rate2f, sizeof(double), N_G, f); fread(sum_rate2b, sizeof(double), N_G, f);
- fread(sum_rate3f, sizeof(double), N_G, f); fread(sum_rate3b, sizeof(double), N_G, f);
- fread(sum_rate4f, sizeof(double), N_G, f); fread(sum_rate4b, sizeof(double), N_G, f);
- fread(sum_rate12f, sizeof(double), N_G, f); fread(sum_rate12b, sizeof(double), N_G, f);
- fread(sum_rate13f, sizeof(double), N_G, f); fread(sum_rate13b, sizeof(double), N_G, f);
- fread(sum_rate17f, sizeof(double), N_G, f); fread(sum_rate17b, sizeof(double), N_G, f);
- fread(sum_rate19f, sizeof(double), N_G, f); fread(sum_rate19b, sizeof(double), N_G, f);
- fread(avg_rate1f, sizeof(double), N_G, f); fread(avg_rate1b, sizeof(double), N_G, f);
- fread(avg_rate2f, sizeof(double), N_G, f); fread(avg_rate2b, sizeof(double), N_G, f);
- fread(avg_rate3f, sizeof(double), N_G, f); fread(avg_rate3b, sizeof(double), N_G, f);
- fread(avg_rate4f, sizeof(double), N_G, f); fread(avg_rate4b, sizeof(double), N_G, f);
- fread(avg_rate12f, sizeof(double), N_G, f); fread(avg_rate12b, sizeof(double), N_G, f);
- fread(avg_rate13f, sizeof(double), N_G, f); fread(avg_rate13b, sizeof(double), N_G, f);
- fread(avg_rate17f, sizeof(double), N_G, f); fread(avg_rate17b, sizeof(double), N_G, f);
- fread(avg_rate19f, sizeof(double), N_G, f); fread(avg_rate19b, sizeof(double), N_G, f);
- fread(sum_electron_density, sizeof(double), N_G, f);
- fread(avg_electron_density, sizeof(double), N_G, f);
- fclose(f);
- printf(">> eduPIC: data loaded : %d electrons %d ions, excited states densities, %d cycles completed before, time is %e [s]\n",N_e,N_i,cycles_done,Time);
- }
- //---------------------------------------------------------------------//
- // save density data //
- //---------------------------------------------------------------------//
- void save_density(void){
- FILE *f;
- double c;
- int m;
- f = fopen("density.dat","w");
- c = 1.0 / (double)(no_of_cycles) / (double)(N_T);
- for(m=0; m<N_G; m++){
- fprintf(f,"%8.5f %12e %12e\n",m * DX, cumul_e_density[m] * c, cumul_i_density[m] * c);
- }
- fclose(f);
- }
- // save Final excited states densities - Artem
- void save_final_excited_densities(void) {
- FILE *f = fopen("excited_densities.dat", "w");
- for(int p=0; p<N_G; p++) {
- fprintf(f, "%8.5f %12e %12e %12e %12e\n", p*DX, e1_density[p], e2_density[p], e3_density[p], e4_density[p]);
- }
- fclose(f);
- }
- //---------------------------------------------------------------------//
- // save EEPF data //
- //---------------------------------------------------------------------//
- void save_eepf(void) {
- FILE *f;
- int i;
- double h,energy;
- h = 0.0;
- for (i=0; i<N_EEPF; i++) {h += eepf[i];}
- h *= DE_EEPF;
- f = fopen("eepf.dat","w");
- for (i=0; i<N_EEPF; i++) {
- energy = (i + 0.5) * DE_EEPF;
- fprintf(f,"%e %e\n", energy, eepf[i] / h / sqrt(energy));
- }
- fclose(f);
- }
- //---------------------------------------------------------------------//
- // save IFED data //
- //---------------------------------------------------------------------//
- void save_ifed(void) {
- FILE *f;
- int i;
- double h_pow,h_gnd,energy;
- h_pow = 0.0;
- h_gnd = 0.0;
- for (i=0; i<N_IFED; i++) {h_pow += ifed_pow[i]; h_gnd += ifed_gnd[i];}
- h_pow *= DE_IFED;
- h_gnd *= DE_IFED;
- mean_i_energy_pow = 0.0;
- mean_i_energy_gnd = 0.0;
- f = fopen("ifed.dat","w");
- for (i=0; i<N_IFED; i++) {
- energy = (i + 0.5) * DE_IFED;
- fprintf(f,"%6.2f %10.6f %10.6f\n", energy, (double)(ifed_pow[i])/h_pow, (double)(ifed_gnd[i])/h_gnd);
- mean_i_energy_pow += energy * (double)(ifed_pow[i]) / h_pow;
- mean_i_energy_gnd += energy * (double)(ifed_gnd[i]) / h_gnd;
- }
- fclose(f);
- }
- //--------------------------------------------------------------------//
- // save XT data //
- //--------------------------------------------------------------------//
- void save_xt_1(xt_distr distr, char *fname) {
- FILE *f;
- int i, j;
- f = fopen(fname,"w");
- for (i=0; i<N_G; i++){
- for (j=0; j<N_XT; j++){
- fprintf(f,"%e ", distr[i][j]);
- }
- fprintf(f,"\n");
- }
- fclose(f);
- }
- void norm_all_xt(void){
- double f1, f2;
- int i, j;
- // normalize all XT data
- f1 = (double)(N_XT) / (double)(no_of_cycles * N_T);
- f2 = WEIGHT / (ELECTRODE_AREA * DX) / (no_of_cycles * (PERIOD / (double)(N_XT)));
- for (i=0; i<N_G; i++){
- for (j=0; j<N_XT; j++){
- pot_xt[i][j] *= f1;
- efield_xt[i][j] *= f1;
- ne_xt[i][j] *= f1;
- ni_xt[i][j] *= f1;
- e1_xt[i][j] *= f1; // Artem
- e2_xt[i][j] *= f1; // Artem
- if (counter_e_xt[i][j] > 0) {
- ue_xt[i][j] = ue_xt[i][j] / counter_e_xt[i][j];
- je_xt[i][j] = -ue_xt[i][j] * ne_xt[i][j] * E_CHARGE;
- meanee_xt[i][j] = meanee_xt[i][j] / counter_e_xt[i][j];
- ioniz_rate_xt[i][j] *= f2;
- } else {
- ue_xt[i][j] = 0.0;
- je_xt[i][j] = 0.0;
- meanee_xt[i][j] = 0.0;
- ioniz_rate_xt[i][j] = 0.0;
- }
- if (counter_i_xt[i][j] > 0) {
- ui_xt[i][j] = ui_xt[i][j] / counter_i_xt[i][j];
- ji_xt[i][j] = ui_xt[i][j] * ni_xt[i][j] * E_CHARGE;
- meanei_xt[i][j] = meanei_xt[i][j] / counter_i_xt[i][j];
- } else {
- ui_xt[i][j] = 0.0;
- ji_xt[i][j] = 0.0;
- meanei_xt[i][j] = 0.0;
- }
- powere_xt[i][j] = je_xt[i][j] * efield_xt[i][j];
- poweri_xt[i][j] = ji_xt[i][j] * efield_xt[i][j];
- }
- }
- }
- void save_all_xt(void){
- char fname[80];
- strcpy(fname,"pot_xt.dat"); save_xt_1(pot_xt, fname);
- strcpy(fname,"efield_xt.dat"); save_xt_1(efield_xt, fname);
- strcpy(fname,"ne_xt.dat"); save_xt_1(ne_xt, fname);
- strcpy(fname,"ni_xt.dat"); save_xt_1(ni_xt, fname);
- strcpy(fname,"je_xt.dat"); save_xt_1(je_xt, fname);
- strcpy(fname,"ji_xt.dat"); save_xt_1(ji_xt, fname);
- strcpy(fname,"powere_xt.dat"); save_xt_1(powere_xt, fname);
- strcpy(fname,"poweri_xt.dat"); save_xt_1(poweri_xt, fname);
- strcpy(fname,"meanee_xt.dat"); save_xt_1(meanee_xt, fname);
- strcpy(fname,"meanei_xt.dat"); save_xt_1(meanei_xt, fname);
- strcpy(fname,"ioniz_xt.dat"); save_xt_1(ioniz_rate_xt, fname);
- strcpy(fname,"e1_xt.dat"); save_xt_1(e1_xt, fname);
- strcpy(fname,"e2_xt.dat"); save_xt_1(e2_xt, fname);
- strcpy(fname,"e3_xt.dat"); save_xt_1(e3_xt, fname);
- strcpy(fname,"e4_xt.dat"); save_xt_1(e4_xt, fname);
- //add
- }
- //---------------------------------------------------------------------//
- // simulation report including stability and accuracy conditions //
- //---------------------------------------------------------------------//
- void check_and_save_info(void){
- FILE *f;
- double plas_freq, meane, kT, debye_length, density, ecoll_freq, icoll_freq, sim_time, e_max, v_max, power_e, power_i, c;
- int i,j;
- bool conditions_OK;
- density = cumul_e_density[N_G / 2] / (double)(no_of_cycles) / (double)(N_T); // e density @ center
- plas_freq = E_CHARGE * sqrt(density / EPSILON0 / E_MASS); // e plasma frequency @ center
- meane = mean_energy_accu_center / (double)(mean_energy_counter_center); // e mean energy @ center
- kT = 2.0 * meane * EV_TO_J / 3.0; // k T_e @ center (approximate)
- sim_time = (double)(no_of_cycles) / FREQUENCY; // simulated time
- ecoll_freq = (double)(N_e_coll) / sim_time / (double)(N_e); // e collision frequency
- icoll_freq = (double)(N_i_coll) / sim_time / (double)(N_i); // ion collision frequency
- debye_length = sqrt(EPSILON0 * kT / density) / E_CHARGE; // e Debye length @ center
- f = fopen("info.txt","w");
- fprintf(f,"########################## eduPIC simulation report ############################\n");
- fprintf(f,"Simulation parameters:\n");
- fprintf(f,"Gap distance = %12.3e [m]\n", L);
- fprintf(f,"# of grid divisions = %12d\n", N_G);
- fprintf(f,"Frequency = %12.3e [Hz]\n", FREQUENCY);
- fprintf(f,"# of time steps / period = %12d\n", N_T);
- fprintf(f,"# of electron / ion time steps = %12d\n", N_SUB);
- fprintf(f,"Voltage amplitude = %12.3e [V]\n", VOLTAGE);
- fprintf(f,"Pressure (Ar) = %12.3e [Pa]\n", PRESSURE);
- fprintf(f,"Temperature = %12.3e [K]\n", T_neutral);
- fprintf(f,"Superparticle weight = %12.3e\n", WEIGHT);
- fprintf(f,"# of simulation cycles in this run = %12d\n", no_of_cycles);
- fprintf(f,"--------------------------------------------------------------------------------\n");
- fprintf(f,"Plasma characteristics:\n");
- fprintf(f,"Electron density @ center = %12.3e [m^{-3}]\n", density);
- fprintf(f,"Plasma frequency @ center = %12.3e [rad/s]\n", plas_freq);
- fprintf(f,"Debye length @ center = %12.3e [m]\n", debye_length);
- fprintf(f,"Electron collision frequency = %12.3e [1/s]\n", ecoll_freq);
- fprintf(f,"Ion collision frequency = %12.3e [1/s]\n", icoll_freq);
- fprintf(f,"--------------------------------------------------------------------------------\n");
- fprintf(f,"Stability and accuracy conditions:\n");
- conditions_OK = true;
- c = plas_freq * DT_E;
- fprintf(f,"Plasma frequency @ center * DT_E = %12.3f (OK if less than 0.20)\n", c);
- if (c > 0.2) {conditions_OK = false;}
- c = DX / debye_length;
- fprintf(f,"DX / Debye length @ center = %12.3f (OK if less than 1.00)\n", c);
- if (c > 1.0) {conditions_OK = false;}
- c = max_electron_coll_freq() * DT_E;
- fprintf(f,"Max. electron coll. frequency * DT_E = %12.3f (OK if less than 0.05)\n", c);
- if (c > 0.05) {conditions_OK = false;}
- c = max_ion_coll_freq() * DT_I;
- fprintf(f,"Max. ion coll. frequency * DT_I = %12.3f (OK if less than 0.05)\n", c);
- if (c > 0.05) {conditions_OK = false;}
- if (conditions_OK == false){
- fprintf(f,"--------------------------------------------------------------------------------\n");
- fprintf(f,"** STABILITY AND ACCURACY CONDITION(S) VIOLATED - REFINE SIMULATION SETTINGS! **\n");
- fprintf(f,"--------------------------------------------------------------------------------\n");
- fclose(f);
- printf(">> eduPIC: ERROR: STABILITY AND ACCURACY CONDITION(S) VIOLATED!\n");
- printf(">> eduPIC: for details see 'info.txt' and refine simulation settings!\n");
- }
- else
- {
- // calculate maximum energy for which the Courant-Friedrichs-Levy condition holds:
- v_max = DX / DT_E;
- e_max = 0.5 * E_MASS * v_max * v_max / EV_TO_J;
- fprintf(f,"Max e- energy for CFL condition = %12.3f [eV]\n", e_max);
- fprintf(f,"Check EEPF to ensure that CFL is fulfilled for the majority of the electrons!\n");
- fprintf(f,"--------------------------------------------------------------------------------\n");
- // saving of the following data is done here as some of the further lines need data
- // that are computed / normalized in these functions
- printf(">> eduPIC: saving diagnostics data\n");
- save_density();
- save_final_excited_densities(); // Artem
- save_eepf();
- save_ifed();
- norm_all_xt();
- save_all_xt();
- fprintf(f,"Particle characteristics at the electrodes:\n");
- fprintf(f,"Ion flux at powered electrode = %12.3e [m^{-2} s^{-1}]\n", N_i_abs_pow * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
- fprintf(f,"Ion flux at grounded electrode = %12.3e [m^{-2} s^{-1}]\n", N_i_abs_gnd * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
- fprintf(f,"Mean ion energy at powered electrode = %12.3e [eV]\n", mean_i_energy_pow);
- fprintf(f,"Mean ion energy at grounded electrode = %12.3e [eV]\n", mean_i_energy_gnd);
- fprintf(f,"Electron flux at powered electrode = %12.3e [m^{-2} s^{-1}]\n", N_e_abs_pow * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
- fprintf(f,"Electron flux at grounded electrode = %12.3e [m^{-2} s^{-1}]\n", N_e_abs_gnd * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
- fprintf(f,"--------------------------------------------------------------------------------\n");
- // calculate spatially and temporally averaged power absorption by the electrons and ions
- power_e = 0.0;
- power_i = 0.0;
- for (i=0; i<N_G; i++){
- for (j=0; j<N_XT; j++){
- power_e += powere_xt[i][j];
- power_i += poweri_xt[i][j];
- }
- }
- power_e /= (double)(N_XT * N_G);
- power_i /= (double)(N_XT * N_G);
- fprintf(f,"Absorbed power calculated as <j*E>:\n");
- fprintf(f,"Electron power density (average) = %12.3e [W m^{-3}]\n", power_e);
- fprintf(f,"Ion power density (average) = %12.3e [W m^{-3}]\n", power_i);
- fprintf(f,"Total power density(average) = %12.3e [W m^{-3}]\n", power_e + power_i);
- fprintf(f,"--------------------------------------------------------------------------------\n");
- fclose(f);
- }
- }
- //------------------------------------------------------------------------------------------//
- // main //
- // command line arguments: //
- // [1]: number of cycles (0 for init) //
- // [2]: "m" turns on data collection and saving //
- //------------------------------------------------------------------------------------------//
- int main (int argc, char *argv[]){
- if (argc == 1) {
- printf(">> eduPIC: error = need starting_cycle argument\n");
- return 1;
- } else {
- strcpy(st0,argv[1]);
- arg1 = atol(st0);
- if (argc > 2) {
- if (strcmp (argv[2],"m") == 0){
- measurement_mode = true; // measurements will be done
- } else {
- measurement_mode = false;
- }
- }
- }
- if (measurement_mode) {
- printf(">> eduPIC: measurement mode: on\n");
- } else {
- printf(">> eduPIC: measurement mode: off\n");
- }
- set_electron_cross_sections_ar();
- set_ion_cross_sections_ar();
- calc_total_cross_sections();
- auto exc = init_excited_distr();
- printf("Excited state population initialized!\n");
- std::ofstream file0("rates.dat"); file0 << std::scientific << std::setprecision(6);
- std::ofstream file1("excited_densities.dat"); file1 << std::scientific << std::setprecision(6);
- //test_cross_sections(); return 1;
- datafile = fopen("conv.dat","a");
- if (arg1 == 0) {
- if (FILE *file = fopen("picdata.bin", "r")) { fclose(file);
- printf(">> eduPIC: Warning: Data from previous calculation are detected.\n");
- printf(" To start a new simulation from the beginning, please delete all output files before running ./eduPIC 0\n");
- printf(" To continue the existing calculation, please specify the number of cycles to run, e.g. ./eduPIC 100\n");
- exit(0);
- }
- no_of_cycles = 1;
- cycle = 1; // init cycle
- init(N_INIT); // seed initial electrons & ions
- printf(">> eduPIC: running initializing cycle\n");
- Time = 0;
- do_one_cycle();
- // std::cout << "CHECKING SUM RATES THROUGH CYCLES ---- " << sum_rate1f[250] << "\n";
- print_excitation_densities();
- for (int i = 0; i < N_G; i++){
- double x = i*DX;
- file0 << avg_rate1f[i] << " " << avg_rate1b[i] << " " << avg_rate2f[i] << " " << avg_rate2b[i] << "\n";
- file1 << x << " " << e1_density[i] << " " << e2_density[i]
- << " " << e3_density[i] << " " << e4_density[i]<< "\n";
- }
- cycles_done = 1;
- } else {
- no_of_cycles = arg1; // run number of cycles specified in command line
- load_particle_data(); // read previous configuration from file
- counter = cycles_done % N_avg;
- printf(">> eduPIC: running %d cycle(s)\n",no_of_cycles);
- for (cycle=cycles_done+1;cycle<=cycles_done+no_of_cycles;cycle++) {
- do_one_cycle();
- // std::cout << "CHECKING SUM RATES THROUGH CYCLES ---- " << sum_rate1f[250] << "|| counter: " << counter << "\n";
- }
- cycles_done += no_of_cycles;
- }
- fclose(datafile);
- save_particle_data();
- if (measurement_mode) {
- check_and_save_info();
- }
- printf(">> eduPIC: simulation of %d cycle(s) is completed.\n",no_of_cycles);
- file0.close();
- file1.close();
- }
Add Comment
Please, Sign In to add comment