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             : 
       9             : #include "PKAGeneratorBase.h"
      10             : #include "MagpieUtils.h"
      11             : #include <algorithm>
      12             : 
      13             : InputParameters
      14         324 : PKAGeneratorBase::validParams()
      15             : {
      16         324 :   InputParameters params = DiscreteElementUserObject::validParams();
      17         324 :   params.addClassDescription("PKA generator user object base class.\n Takes pdf and samples PKAs "
      18             :                              "due to various interactions.");
      19         324 :   return params;
      20           0 : }
      21             : 
      22         179 : PKAGeneratorBase::PKAGeneratorBase(const InputParameters & parameters)
      23         179 :   : DiscreteElementUserObject(parameters)
      24             : {
      25         179 :   setRandomResetFrequency(EXEC_TIMESTEP_END);
      26         179 : }
      27             : 
      28             : void
      29      718143 : PKAGeneratorBase::setPosition(MyTRIM_NS::IonBase & ion) const
      30             : {
      31      718143 :   ion._pos = MagpieUtils::randomElementPoint(*_current_elem, getRandomPoint());
      32      718143 : }
      33             : 
      34             : int
      35     1355334 : PKAGeneratorBase::ionTag(const MyTRIMRasterizer::PKAParameters & pka_parameters,
      36             :                          Real Z,
      37             :                          Real m) const
      38             : {
      39             :   // this function relies on the exact representation of whole numbers in IEEE floating point
      40             :   // numbers up to a reasonable upper limit [Z < m < 300]
      41             : 
      42             :   const auto & mZ = pka_parameters._mass_charge_tuple;
      43             : 
      44             :   // element not found in rasterizer table
      45     1355334 :   if (pka_parameters._index_Z[Z].first == 0)
      46             :     return -1;
      47             : 
      48             :   // only one isotope of this element is present
      49       62022 :   if (pka_parameters._index_Z[Z].first == 1)
      50             :   {
      51       62022 :     auto & t = mZ[pka_parameters._index_Z[Z].second];
      52       62022 :     if (std::abs(m - std::get<0>(t)) < std::get<2>(t))
      53       62019 :       return pka_parameters._index_Z[Z].second;
      54             :     else
      55             :       return -1;
      56             :   }
      57             : 
      58             :   // start the search at the first matching Z
      59           0 :   for (auto i = pka_parameters._index_Z[Z].second; i < mZ.size(); ++i)
      60           0 :     if (std::get<1>(mZ[i]) == Z && std::abs(m - std::get<0>(mZ[i])) < std::get<2>(mZ[i]))
      61           0 :       return i;
      62             : 
      63             :   // no matching mass (isotope) found for the given Z
      64             :   return -1;
      65             : }
      66             : 
      67             : void
      68      727173 : PKAGeneratorBase::setRandomDirection(MyTRIM_NS::IonBase & ion) const
      69             : {
      70             :   Real nsq, x1, x2;
      71             : 
      72             :   // Marsaglia's method for uniformly sampling the surface of the sphere
      73             :   do
      74             :   {
      75      924711 :     x1 = 2 * getRandomReal() - 1.0;
      76      924711 :     x2 = 2 * getRandomReal() - 1.0;
      77      924711 :     nsq = x1 * x1 + x2 * x2;
      78      924711 :   } while (nsq >= 1);
      79             : 
      80             :   // construct normalized direction vector
      81      727173 :   ion._dir(0) = 2.0 * x1 * std::sqrt(1.0 - nsq);
      82      727173 :   ion._dir(1) = 2.0 * x2 * std::sqrt(1.0 - nsq);
      83      727173 :   ion._dir(2) = 1.0 - 2.0 * nsq;
      84      727173 : }