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 "PKAGeneratorAlphaDecay.h"
      10             : #include "MagpieUtils.h"
      11             : #include <algorithm>
      12             : 
      13             : registerMooseObject("MagpieApp", PKAGeneratorAlphaDecay);
      14             : 
      15             : InputParameters
      16           9 : PKAGeneratorAlphaDecay::validParams()
      17             : {
      18           9 :   InputParameters params = PKAGeneratorBase::validParams();
      19           9 :   params.addClassDescription("A PKAGenerator for starting alpha particles from decay\nDecay data "
      20             :                              "is retrieved from file data/alpha_decay/alpha_decay.txt");
      21          18 :   params.addRequiredParam<std::vector<unsigned int>>("ZAID",
      22             :                                                      "The Z/A ids corresponding to the var "
      23             :                                                      "arguments in the rasterizer [1e4 * A + 10 * "
      24             :                                                      "Z + state]");
      25          18 :   MooseEnum timeUnit("second=0 millisecond=1 microsecond=2", "second");
      26          18 :   params.addParam<MooseEnum>("time_unit", timeUnit, "Unit of time used in this model");
      27           9 :   return params;
      28           9 : }
      29             : 
      30           5 : PKAGeneratorAlphaDecay::PKAGeneratorAlphaDecay(const InputParameters & parameters)
      31          15 :   : PKAGeneratorBase(parameters), _zaids(getParam<std::vector<unsigned int>>("ZAID"))
      32             : {
      33          10 :   switch (getParam<MooseEnum>("time_unit"))
      34             :   {
      35           0 :     case 0:
      36           0 :       _time_conversion = 1.0;
      37           0 :       break;
      38           0 :     case 1:
      39           0 :       _time_conversion = 1.0e3;
      40           0 :       break;
      41           5 :     case 2:
      42           5 :       _time_conversion = 1.0e6;
      43           5 :       break;
      44             :   }
      45           5 :   readAlphaData();
      46           5 : }
      47             : 
      48             : void
      49           5 : PKAGeneratorAlphaDecay::readAlphaData()
      50             : {
      51           5 :   auto path = std::getenv("ALPHA_DIR");
      52           5 :   if (path == NULL)
      53           0 :     mooseError("Set $ALPHA_DIR to the directory holding the alpha_decay.txt file.");
      54             : 
      55           5 :   std::string s = "alpha_decay.txt";
      56           5 :   std::string filename = path + s;
      57             : 
      58           5 :   if (!MooseUtils::checkFileReadable(filename, false, false))
      59           0 :     mooseError("Missing/Corrupted ", filename, " file");
      60             : 
      61           5 :   std::ifstream infile(filename.c_str());
      62             :   unsigned int nentries;
      63             :   infile >> nentries;
      64        2395 :   for (unsigned int j = 0; j < nentries; ++j)
      65             :   {
      66             :     unsigned int zaid;
      67             :     unsigned int n_modes;
      68             :     infile >> zaid >> n_modes;
      69        2390 :     std::vector<DecayData> decay(n_modes, DecayData());
      70       12760 :     for (unsigned l = 0; l < n_modes; ++l)
      71             :     {
      72             :       unsigned int Z, A;
      73             :       Real half_life, energy, intensity;
      74             :       infile >> Z >> A >> half_life >> energy >> intensity;
      75             :       // adjust half_life unit; file units are
      76       10370 :       decay[l]._decay_constants = std::log(2.0) / (half_life * _time_conversion);
      77       10370 :       decay[l]._alpha_energies = energy;
      78       10370 :       decay[l]._intensities = intensity;
      79             :     }
      80             : 
      81        2390 :     if (std::find(_zaids.begin(), _zaids.end(), zaid) != _zaids.end())
      82             :     {
      83          10 :       if (_decay_data_sets.find(zaid) != _decay_data_sets.end())
      84           0 :         mooseError("Identical Z/A id was provided twice");
      85          10 :       _decay_data_sets[zaid] = decay;
      86             :     }
      87             :   }
      88           5 :   infile.close();
      89             : 
      90             :   // complete the data by dealing with absent nuclides ==> stable
      91          25 :   for (auto & z : _zaids)
      92             :   {
      93          20 :     if (_decay_data_sets.find(z) == _decay_data_sets.end())
      94             :     {
      95          10 :       std::vector<DecayData> decay(1, DecayData());
      96          10 :       decay[0]._decay_constants = 0.0;
      97          10 :       decay[0]._alpha_energies = 0.0;
      98          10 :       decay[0]._intensities = 1.0;
      99          10 :       _decay_data_sets[z] = decay;
     100             :     }
     101             :   }
     102             : 
     103             : // in debug mode we want some info on what is stored in _decay_data_sets
     104             : #if DEBUG
     105             :   _console << "\nAlpha decay read from file\n";
     106             :   for (auto & d : _decay_data_sets)
     107             :   {
     108             :     _console << "ZAID: " << d.first << " has " << d.second.size() << " decay channels:\n";
     109             :     for (unsigned int j = 0; j < d.second.size(); ++j)
     110             :       _console << "lambda = " << d.second[j]._decay_constants << " "
     111             :                << "energy = " << d.second[j]._alpha_energies << " "
     112             :                << "intensity = " << d.second[j]._intensities << "\n";
     113             :   }
     114             :   _console << std::endl;
     115             : #endif
     116          10 : }
     117             : 
     118             : void
     119       24000 : PKAGeneratorAlphaDecay::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_list,
     120             :                                    const MyTRIMRasterizer::PKAParameters & pka_parameters,
     121             :                                    const MyTRIMRasterizer::AveragedData & averaged_data) const
     122             : {
     123       24000 :   const auto dt = pka_parameters._dt;
     124       24000 :   const auto vol = pka_parameters._volume;
     125       24000 :   const auto recoil_rate_scaling = pka_parameters._recoil_rate_scaling;
     126             : 
     127             :   mooseAssert(dt >= 0,
     128             :               "Passed a negative time window into PKAFissionFragmentNeutronics::appendPKAs");
     129             :   mooseAssert(vol >= 0, "Passed a negative volume into PKAFissionFragmentNeutronics::appendPKAs");
     130             : 
     131       24000 :   if (averaged_data._elements.size() != _zaids.size())
     132           0 :     mooseError("Size of averaged_data and ZAID must be equal");
     133             : 
     134      120000 :   for (unsigned int nuclide = 0; nuclide < _zaids.size(); ++nuclide)
     135             :   {
     136       96000 :     unsigned int zaid = _zaids[nuclide];
     137             : 
     138             :     auto it = _decay_data_sets.find(zaid);
     139             : 
     140       96000 :     if (it == _decay_data_sets.end())
     141           0 :       mooseError("ZAID ", zaid, " not found in decay data set");
     142             : 
     143             :     auto & decay = it->second;
     144      312000 :     for (unsigned int l = 0; l < decay.size(); ++l)
     145             :     {
     146      216000 :       unsigned int num_decay = std::floor(recoil_rate_scaling * vol * decay[l]._intensities *
     147      216000 :                                               averaged_data._elements[nuclide] *
     148      216000 :                                               (1.0 - std::exp(-decay[l]._decay_constants * dt)) +
     149      216000 :                                           getRandomReal());
     150             : 
     151      216174 :       for (unsigned i = 0; i < num_decay; ++i)
     152             :       {
     153         174 :         std::vector<MyTRIM_NS::IonBase> ion(2);
     154             : 
     155             :         // set stopping criteria
     156         174 :         ion[0].setEf();
     157         174 :         ion[1].setEf();
     158             : 
     159             :         // get parent Z/A
     160         174 :         unsigned int parent_Z = MagpieUtils::getZFromZAID(zaid);
     161         174 :         unsigned int parent_A = MagpieUtils::getAFromZAID(zaid);
     162             : 
     163             :         // set Z/A values; [0] is alpha, [1] is recoil
     164         174 :         ion[0]._Z = 2;
     165         174 :         ion[0]._m = 4.0;
     166             : 
     167         174 :         ion[1]._Z = parent_Z - 2;
     168         174 :         ion[1]._m = parent_A - 4;
     169             : 
     170         174 :         ion[0]._tag = ionTag(pka_parameters, ion[0]._Z, ion[0]._m);
     171         174 :         ion[1]._tag = ionTag(pka_parameters, ion[1]._Z, ion[1]._m);
     172             : 
     173             :         // set the energies; use momentum conservation
     174         174 :         ion[0]._E = decay[l]._alpha_energies;
     175         174 :         ion[1]._E = decay[l]._alpha_energies * ion[0]._m / ion[1]._m;
     176             : 
     177             :         // set location of the fission event
     178         174 :         setPosition(ion[0]);
     179         174 :         ion[1]._pos = ion[0]._pos;
     180             : 
     181             :         // set random direction for ion 1 and opposite direction for ion 2
     182         174 :         setRandomDirection(ion[0]);
     183         174 :         ion[1]._dir = -ion[0]._dir;
     184             : 
     185             :         // add PKAs to list
     186         174 :         ion_list.push_back(ion[0]);
     187         174 :         ion_list.push_back(ion[1]);
     188         174 :       }
     189             :     }
     190             :   }
     191       24000 : }