Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : #ifdef GSL_ENABLED
       9             : 
      10             : #pragma once
      11             : 
      12             : #include <gsl/gsl_errno.h>
      13             : #include <gsl/gsl_matrix.h>
      14             : #include <gsl/gsl_odeiv2.h>
      15             : 
      16             : #include <mytrim/material.h>
      17             : 
      18             : class SimconfType;
      19             : class IonBase;
      20             : 
      21             : enum nrt_type
      22             : {
      23             :   ENERGY,
      24             :   TOTAL,
      25             :   NET,
      26             :   NET_DERIVATIVE
      27             : };
      28             : 
      29             : /**
      30             :  * Implements the computation of polyatomic displacement functions using Parkin & Coulter's
      31             :  * method
      32             :  */
      33             : class PolyatomicDisplacementFunctionBase
      34             : {
      35             : public:
      36             :   /// default constructor
      37             :   PolyatomicDisplacementFunctionBase(std::vector<MyTRIM_NS::Element> polyatomic_material,
      38             :                                      nrt_type damage_function_type,
      39             :                                      std::vector<std::vector<Real>> Ecap);
      40             : 
      41             :   virtual ~PolyatomicDisplacementFunctionBase();
      42             : 
      43             :   static int gslInterface(Real energy, const Real disp[], Real f[], void * params);
      44             : 
      45             :   virtual std::vector<Real> getRHS(Real energy) = 0;
      46             : 
      47             :   /// computes the displacement function from current last energy to Emax
      48             :   void advanceDisplacements(Real Emax);
      49             : 
      50             :   /// computes the integral of the displacement function
      51             :   void computeDisplacementFunctionIntegral();
      52             : 
      53             :   /// obtain quadrature rule from GSL
      54             :   static void gslQuadRule(unsigned int quad_order,
      55             :                           std::vector<Real> & quad_points,
      56             :                           std::vector<Real> & quad_weights);
      57             : 
      58             :   ///@{ some getters needed for accessing this pointer in odeRHS
      59     2340391 :   unsigned int nSpecies() const { return _n_species; }
      60             :   unsigned int problemSize() const { return _problem_size; }
      61       47538 :   unsigned int nEnergySteps() const { return _energy_history.size(); }
      62             :   unsigned int findSpeciesIndex(int atomic_number, Real mass_number) const;
      63          42 :   Real minEnergy() const { return _energy_history[0]; }
      64      140788 :   Real energyPoint(unsigned int na) const { return _energy_history[na]; }
      65       13440 :   Real lambda(unsigned int i, unsigned int j) const { return _lambda[i][j]; }
      66      784765 :   Real numberFraction(unsigned int i) const { return _material->_element[i]._t; }
      67             :   Real numberDensity(unsigned int i) const { return _material->_element[i]._t * _material->_arho; }
      68             :   ///@}
      69             : 
      70             :   ///@{ linear interpolation of the damage function
      71             :   Real
      72             :   linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const;
      73             :   Real linearInterpolationFromFlatIndex(Real energy, unsigned int index) const;
      74             :   ///@}
      75             : 
      76             :   /// linear interpolation of the integral of the damage function
      77             :   Real linearInterpolationIntegralDamageFunction(Real energy,
      78             :                                                  unsigned int i,
      79             :                                                  unsigned int j = 0,
      80             :                                                  unsigned int l = 0) const;
      81             : 
      82             :   /// gets stopping power for a given species and energy; non-const because it uses _ions so no need to construct ion
      83             :   Real stoppingPower(unsigned int species, Real energy);
      84             : 
      85             :   /// returns derivative of the stopping power of species w.r.t. to changes in the number fraction of changing_species
      86             :   Real stoppingPowerDerivative(unsigned int species, unsigned int changing_species, Real energy);
      87             : 
      88             :   /// retrives the grid index for an energy value
      89             :   unsigned int energyIndex(Real energy) const;
      90             : 
      91             : protected:
      92             :   /// override the mapIndex function: flattens ijl to single index n
      93             :   unsigned int mapIndex(unsigned int i, unsigned int j = 0, unsigned int l = 0) const;
      94             : 
      95             :   /// computes the integral int_0^t dT T * d(sigma_ij) / dT for species combination i, j and small t
      96             :   Real
      97             :   weightedScatteringIntegral(Real energy, Real energy_limit, unsigned int i, unsigned int j) const;
      98             : 
      99             :   /// atomic scattering cross section d sigma_ij / dT, i: projectile, j: target
     100             :   Real
     101             :   scatteringCrossSection(unsigned int i, unsigned int j, Real energy, Real recoil_energy) const;
     102             : 
     103             :   /// the universal scattering function by Lindhard usually denoted as f(t^1/2)
     104             :   Real universalF(Real xi) const;
     105             : 
     106             :   /// displacement probability
     107             :   virtual Real displacementProbability(unsigned int k, Real energy) const;
     108             : 
     109             :   /// capture probability
     110             :   virtual Real
     111             :   nonCaptureProbability(unsigned int i, unsigned int k, Real energy, Real recoil_energy) const;
     112             : 
     113             :   /// a helper function called from linearInterpolation
     114             :   Real linearInterpolationHelper(Real energy, unsigned int energy_index, unsigned int index) const;
     115             : 
     116             :   /// damage function type [nij and gij, respectively in PK JNM 101, 1981; or nu_i JNM 88, (1980)]
     117             :   nrt_type _damage_function_type;
     118             : 
     119             :   /// order of the quadrature
     120             :   unsigned int _quad_order;
     121             : 
     122             :   /// the number of different species in the material
     123             :   unsigned int _n_species;
     124             : 
     125             :   /// the number of species indices of this object
     126             :   unsigned int _n_indices = 0;
     127             : 
     128             :   /// the size of the problem = _n_species**2 for TOTAL & NET, _n_species for ENERGY
     129             :   unsigned int _problem_size;
     130             : 
     131             :   /// the current maximum energy to which PC equations are solved
     132             :   std::vector<Real> _energy_history;
     133             : 
     134             :   /// The displacement values nu_ij flattened into 1D array
     135             :   std::vector<std::vector<Real>> _displacement_function;
     136             : 
     137             :   /// The integral of the displacement function over energy
     138             :   std::vector<std::vector<Real>> _displacement_function_integral;
     139             : 
     140             :   /// Ecap_ik average residual energy of type i atom to not be trapped at k-site
     141             :   std::vector<std::vector<Real>> _Ecap;
     142             : 
     143             :   /// _lambda = 4.0 * Ai * Aj / (Ai + Aj) / (Ai + Aj)
     144             :   std::vector<std::vector<Real>> _lambda;
     145             : 
     146             :   /// Bohr radius in A
     147             :   const Real _abohr = 0.529177;
     148             : 
     149             :   /**
     150             :    * energy threshold below which the integral over cross sections is treated using
     151             :    * the asymptotic form of the cross sections for low energies and nu_j(T) is replaced
     152             :    * by nu_j(T) \approx T.
     153             :    */
     154             :   Real _asymptotic_threshold = 0.1;
     155             : 
     156             :   /// integral of the scattering cross sections int_{Edisp}^{Lambda_ij E} dT d(sigma_ij)/dT; Edisp: target displacement threshold
     157             :   // std::vector<std::vector<Real>> _scattering_xs_integral;
     158             : 
     159             :   ///@{ MYTRIM objects to use compute stopping power and cross sections
     160             :   std::unique_ptr<MyTRIM_NS::SimconfType> _simconf;
     161             :   std::unique_ptr<MyTRIM_NS::MaterialBase> _material;
     162             :   std::vector<std::unique_ptr<MyTRIM_NS::IonBase>> _ions;
     163             :   ///@}
     164             : 
     165             :   ///@{ GSL ODE system and driver to solve ODE
     166             :   gsl_odeiv2_system _sys;
     167             :   gsl_odeiv2_driver * _ode_driver;
     168             :   ///@}
     169             : 
     170             :   ///@{ Gaussian quadrature rule
     171             :   std::vector<Real> _quad_points;
     172             :   std::vector<Real> _quad_weights;
     173             :   ///@}
     174             : };
     175             : 
     176             : #endif // GSL_ENABLED