Split eii

.gitignore

@@ -43,7 +43,7 @@ core
 scripts/scattering/R_plots
 *.pickle
 scripts/scattering/dmg_data
-valgrind-out*
+valgrind*
 input/_batches/*
 scripts/scattering/reflections/*
 scripts/scattering/scalepack/*

include/ComputeRateParam.h

@@ -63,7 +63,7 @@ class ComputeRateParam
 	~ComputeRateParam();
 
 	// Halfwidth = 5/Constant::Time -> 5 fs half width.
-	RateData::Atom SolveAtomicRatesAndPlasmaBEB(vector<int> Max_occ, vector<int> Final_occ, vector<bool> shell_check, bool calculate_secondary_ionisation, Grid &Lattice, vector<RadialWF> &Orbitals, Potential &U, Input & Inp, ofstream & log);
+	RateData::Atom SolveAtomicRatesAndPlasmaBEB(vector<int> Max_occ, vector<int> Final_occ, vector<bool> shell_check, bool calculate_secondary_ionisation, ofstream & log);
 	// // Atomic.
 	// int SetupAndSolve(ofstream & log);
 	// // Molecular.

include/Constant.h

@@ -194,7 +194,8 @@ namespace RateData {
 
 	bool ReadRates(const string & input, vector<RateData::Rate> & PutHere);
 	bool ReadEIIParams(const string & input, vector<RateData::EIIdata> & PutHere);
-	bool InterpolateRates(const string & rate_location, const string & rate_file_type,vector<RateData::Rate> & PutHere, double photon_energy,double allowed_interp = 1000);
+	bool ReadDecayRates(const string & rate_location, const string & rate_file_type, vector<RateData::Rate> & PutHere,int num_allowed_configs);
+	bool InterpolateRates(const string & rate_location, const string & rate_file_type, vector<RateData::Rate> & PutHere, double photon_energy,double allowed_interp = 1000);
 	void WriteRates(const string& fname, const vector<RateData::Rate>& rateVector);
 	void WriteEIIParams(const string& fname, const vector<RateData::EIIdata>& eiiVector);
 }

include/FreeDistribution.h

@@ -49,28 +49,47 @@ struct indexed_knot
  * @brief Electron distribution class.
  * @details Represents a statistical distribution of electron density. Internal units are atomic units.
  * @note F is referred to as the distribution throughout the code. F[i] returns the i'th spline factor, but
- * F.f is the container that holds the spline factors. F(e) expands out the basis to return the density at energy e divided by local knot width.
+ * F.f_array[0], (in prior versions F.f), is the container that holds the spline factors for the total continuum (later indices are for tracking cascades initiated by each atomic species). F(e) expands out the basis to return the density at energy e divided by local knot width.
  */
 class Distribution
 {
 public:
     Distribution() {
-        f.resize(size);   
+        f_array.resize(num_continuums);
+        for (auto& f : f_array){
+            f.resize(size);}   
     }
 
-    /**
-     * @brief Get spline factor at 
-     * @param n 
-     * @return 
-     */
-    inline double& operator[](size_t n) {
-        return this->f[n];
+
+    inline std::vector<double>& operator[](size_t _c) {
+        return this->f_array[_c];
     }
 
-    inline double operator[](size_t n) const{
-        return this->f[n];
+    inline std::vector<double> operator[](size_t _c) const{
+        return this->f_array[_c];
     }
 
+    // /**
+    //  * @brief Get spline factors at 
+    //  * @param n 
+    //  * @return 
+    //  */
+    // inline double& operator()(size_t _c, size_t n) {
+    //     return this->f_array[_c][n];
+    // }
+    // inline double operator()(size_t _c, size_t n) const{
+    //     return this->f_array[_c][n];
+    // }
+
+    // // Accesses the continuum, not atom-split parts of it.
+    // inline double& operator[](size_t n) {
+    //     return this->f_array[0][n];
+    // }
+
+    // inline double operator[](size_t n) const{
+    //     return this->f_array[0][n];
+    // }
+
     // vector-space algebra
 
     /**
@@ -79,8 +98,10 @@ class Distribution
      * @return Distribution& 
      */
     Distribution& operator+=(const Distribution& d) {
-        for (size_t i=0; i<size; i++) {
-            f[i] += d.f[i];
+        for (size_t _c = 0; _c < num_continuums; _c++){  // Iterate through each continuum
+            for (size_t i=0; i<size; i++) {
+                f_array[_c][i] += d.f_array[_c][i];
+            }
         }
         // total += d.total;
         return *this;
@@ -92,8 +113,10 @@ class Distribution
      * @return Distribution& 
      */
     Distribution& operator*=(double x) {
-        for (size_t i=0; i<size; i++) {
-            f[i] *= x;
+        for (auto& f : this->f_array){
+            for (size_t i=0; i<size; i++) {
+                f[i] *= x;
+            }
         }
         // total *= d.total;
         return *this;
@@ -105,7 +128,10 @@ class Distribution
      * @return Distribution& 
      */
     Distribution& operator=(const Distribution& d) {
-        f = d.f;
+        f_array = d.f_array;
+        // for (size_t _c = 0; _c < num_continuums; _c++){ 
+        //     f_array[_c] = d.f_array[_c];
+        // }
         // total = d.total;
         return *this;
     }
@@ -117,9 +143,11 @@ class Distribution
      */
     Distribution& operator=(double y) {
         // total = y;
-        f.resize(size);
-        for (size_t i=0; i<size; i++) {
-            f[i] = y;
+        for (auto& f : this->f_array){
+            f.resize(size);
+            for (size_t i=0; i<size; i++) {
+                f[i] = y;
+            }
         }
         return *this;
     }
@@ -146,13 +174,13 @@ class Distribution
      * @return Distribution& 
      * @todo make clearer this is only for static grid case or just put into one function.
      */
-    void set_distribution_STATIC_ONLY(vector<double> new_knot, vector<double> new_f);
-    void set_spline_factors(vector<double> new_f);
+    void set_distribution_STATIC_ONLY(size_t _c, vector<double> new_knot, vector<double> new_f);
+    void set_spline_factors(size_t _c, vector<double> new_f);
     // loads the knot and appropriately updates params like the overlap matrix as done in set_parameters.
     static void load_knot(vector<double> loaded_knot);
 
-    double norm() const;
-    double integral() const; //unused
+    double norm(size_t _c) const;
+    //double integral(size_t _c) const; //unused
     // modifiers
 
     /**
@@ -162,36 +190,36 @@ class Distribution
      * This seems correct, given that dfdt is inputted for the spline basis's Sinv (S_inverse * <VectorXd input>) function, which names the input deltaf.
      * @param v deltaf
      */
-    void applyDeltaF(const Eigen::VectorXd& dfdt, const int & threads);
-
+    void applyDeltaF(size_t a, const Eigen::VectorXd& dfdt, const int & threads);
+    void applyDeltaF_element_scaled(size_t a, const Eigen::VectorXd& dfdt, const int & threads);
 
     // Q functions
     // Computes the dfdt vector v based on internal f
     // e.g. dfdt v; F.calc_Qee(v);
-    void get_Q_eii (Eigen::VectorXd& v, size_t a, const bound_t& P, const int & threads) const;
-    void get_Q_tbr (Eigen::VectorXd& v, size_t a, const bound_t& P, const int & threads) const;
-    void get_Q_ee  (Eigen::VectorXd& v, const int & threads) const;
+    void get_Q_eii (size_t _c, Eigen::VectorXd& v, size_t a, const bound_t& P, const int & threads) const;
+    void get_Q_tbr (size_t _c, Eigen::VectorXd& v, size_t a, const bound_t& P, const int & threads) const;
+    void get_Q_ee  (size_t _c, Eigen::VectorXd& v, const int & threads) const;
 
     void get_Jac_ee (Eigen::MatrixXd& J) const; // Returns the Jacobian of Qee
 
     /// N is the Number density (inverse au^3) of particles to be added at energy e.
     static void addDeltaLike(Eigen::VectorXd& v, double e, double N);
     /// Adds a Dirac delta to the distribution
-    void addDeltaSpike(double N, double e);
+    void addDeltaSpike(size_t a, double N, double e);
     /// Applies the loss term to the distribution 
-    void addLoss(const Distribution& d, const LossGeometry& l, double charge_density);
-    void addFiltration(const Distribution& d, const Distribution& bg,const LossGeometry &l);
+    void addLoss(size_t _c, const Distribution& d, const LossGeometry& l, double charge_density);
+    void addFiltration(size_t _c, const Distribution& d, const Distribution& bg,const LossGeometry &l);
 
     /// Sets the object to have a MB distribution
-    void add_maxwellian(double N, double T);
+    void add_maxwellian(size_t _c, double N, double T);
 
     /**
      * @brief 
      * 
      * @param densities Densities corresponding to each energy.
      * @param energies Energies in ascending order.
      */
-    void add_density_distribution(std::vector<vector<double>>);
+    void add_density_distribution(size_t _c, std::vector<vector<double>>);
 
     // Precalculators
     static void Gamma_eii( eiiGraph& Gamma, const std::vector<RateData::EIIdata>& eii, size_t J) {
@@ -228,9 +256,9 @@ class Distribution
      * @param num_pts 
      * @return 
      */
-    std::string output_densities(size_t num_pts, std::vector<double> reference_knots) const;
+    std::string output_densities(size_t _c, size_t num_pts, std::vector<double> reference_knots) const;
 
-    std::vector<double> get_densities(size_t num_pts, std::vector<double> reference_knots) const;
+    std::vector<double> get_densities(size_t _c, size_t num_pts, std::vector<double> reference_knots) const;
     static std::string output_knots_eV();
     // 
     /**
@@ -244,7 +272,7 @@ class Distribution
     // (Unused) This does electron-electron because it is CURSED
     void from_backwards_Euler(double dt, const Distribution& prev_step, double tolerance, unsigned maxiter);
 
-    double operator()(double e) const;
+    double operator()(size_t _c, double e) const;
 
     /**
      * @brief The setup function.
@@ -277,20 +305,20 @@ class Distribution
         return basis.mb_max_over_kT;
     }
     static size_t size;
-    double my_size(){return f.size();}
+    static size_t num_continuums;
     static std::vector<double> load_knots_from_history(size_t step_idx);
     static size_t most_recent_knot_change_idx(size_t step_idx);
     static size_t next_knot_change_idx(size_t step_idx);
     static std::vector<double> get_knots_from_history(size_t step_idx);
     static void set_knot_history(size_t i, std::vector<double> replacement_knot){knots_history[i]={i,replacement_knot};}
     // history of grid points (for dynamic grid)
     static std::vector<indexed_knot> knots_history;  
-    int container_size(){return f.size();}  
+    int container_size(){return f_array[0].size();}  
 
     static bool reset_on_next_grid_update; //TODO part of a duct tape implementation of FIND_INITIAL_DIRAC
 
 private:
-    std::vector<double> f;  // Spline expansion factors
+    std::vector<std::vector<double>> f_array;  // Spline expansion factors, placed in vector so have option to track a continuum for each element's cascades.
     static SplineIntegral basis;
     static size_t CoulombLog_cutoff;
     /// Coulomb repulsion is ignored if density is below this threshold

include/HybridIntegrator.hpp

@@ -351,7 +351,7 @@ void Hybrid<T>::step_stiff_part(unsigned n){
             y_transient[next_rel_idx] = tmp;
             y_transient[next_rel_idx] += dydt;
             prev += y_transient[next_rel_idx];
-            diff = prev.norm()/y_transient[next_rel_idx].norm(); // Seeking convergence
+            diff = prev.norm(0)/y_transient[next_rel_idx].norm(0); // Seeking convergence
             idx++;
         }
         if(idx==stiff_max_iter){
@@ -439,11 +439,11 @@ void Hybrid<T>::initialise_transient_y(int n) {  // n is the last calculated ste
         }
     }
     // sample to catch the error of y(n) != y_transient(mini_n).
-    assert(this->y[n].F[0] == y_transient[mini_n].F[0]);
+    assert(this->y[n].F[0][0] == y_transient[mini_n].F[0][0]);
     if (this->y[n].F.container_size() > 3)
-        assert(this->y[n].F[3] == y_transient[mini_n].F[3]);
+        assert(this->y[n].F[0][3] == y_transient[mini_n].F[0][3]);
     if (this->y[n].F.container_size() > 5)
-        assert(this->y[n].F[5] == y_transient[mini_n].F[5]);
+        assert(this->y[n].F[0][5] == y_transient[mini_n].F[0][5]);
 }
 
 
@@ -566,7 +566,7 @@ void Hybrid<T>::initialise_transient_y_v2(int n){
             y_transient[mini_n-i] = this->y[n-i];
         }
     }
-    assert(this->y[n].F[3] == y_transient[(mini_n)%(this->order)].F[3]);  
+    assert(this->y[n].F[0][3] == y_transient[(mini_n)%(this->order)].F[0][3]);  
 }
 
 
@@ -596,7 +596,7 @@ void Hybrid<T>::backward_Euler(unsigned n){
 
         tmp *= -1;
         tmp += this->y[n+1];
-        diff = tmp.norm()/this->y[n].norm();
+        diff = tmp.norm(0)/this->y[n].norm(0);
         idx++;
     }
     this->y[n+1] += old;

include/RateSystem.h

@@ -60,7 +60,7 @@ class state_type
     state_type& operator=(const double x);
     // state_type& operator=(const state_type& s2);
 
-    double norm() const;
+    double norm(size_t _c) const;
 
     // Defines number and style of atomP
     // Resizes the container to fit all of the states present in the atom ensemble

include/config.h

@@ -22,7 +22,7 @@ const int GLOBAL_BSPLINE_ORDER = 3;  // 1 = rectangles, 2=linear, 3=quadratic  A
 /// Disable features
 //#define NO_PLOTTING // Turns off live saves of the free-electron energy distribution to _live_plot.png. Disables use of python 
 //#define NO_BACKUP_SAVING // Disables the hourly saves of the data to  
-
+//#define TRACK_SINGLE_CONTINUUM // Just track the total electron density. Computationally expensive to turn off. If off, tracks one electron distribution for each species (element) defined in input file, tracking the primary electrons that they release and the secondary electrons those electrons free from ALL species.
 
 /// Asynchronous solver
 //#define NO_MINISTEPS   // Disables the asynchronous implementation of the solver, stepping the free (E-E) and bound (everything else) solvers together.

input/I3C/I3C_15kev.mol

@@ -0,0 +1,46 @@
+#ATOMS
+C 1000
+N 125
+O 500
+I_fast 375
+
+#BOUND_FREE_EXCLUSIONS
+I_fast
+
+#VOLUME
+40000        // Volume per molecule in Angstrom^3.
+500         // Length scale of sample in Angstrom. Used for effective escape rate of photo-electrons (assumes surrounded by a void).
+none // Spatial boundary shape - options are none, spherical, cylindrical, planar
+
+#PULSE
+15000         // Photon energy in eV.
+10           // Pulse width in femtoseconds (defined as FWHM for Gaussian pulse).
+gaussian       // Pulse shape
+7.2  // total sim time in fwhm
+1.2 // time before pulse peak to start at in fwhm 
+
+#USE_COUNT
+true    // active (using count)?
+0.1  // Photon density (x10^12 ph.µm-2) (3.5e12 um-2)
+
+#NUMERICAL
+2000      // Initial guess for number of time step points for rate equation.
+26            // Number of threads in OpenMP.
+
+#DYNAMIC_GRID
+M              // Grid regions preset, options are: dismal, low, medium, high, (referring to accuracy),  and no_dirac (static region is used rather than dynamic dirac regions)
+5.05            // Grid update period in fs, (dynamic grid only).
+
+#OUTPUT
+800         // Number of time steps in the output files.
+4000        // number of free-electron grid points in output file.
+N            // Write atomic charges in a separate file (Y/N)?
+Y            // Write intensity in a separate file (Y/N)?
+Y            // Write data for molecular dynamics (MD) in a separate file (Y/N)?
+
+#DEBUG
+99 //10     // Simulation cutoff time in fs
+0.01    // Interval to update current timestep [fs].
+####END#####
+
+Notes: