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 "DiscreteFissionPKAPDF.h"
      10             : #include "MooseError.h"
      11             : #include "MooseRandom.h"
      12             : #include "MagpieUtils.h"
      13             : #include "MooseUtils.h"
      14             : 
      15             : #include "libmesh/utility.h"
      16             : 
      17             : #include <fstream>
      18             : #include <map>
      19             : 
      20           0 : DiscreteFissionPKAPDF::DiscreteFissionPKAPDF()
      21             :   : DiscretePKAPDFBase(),
      22           0 :     _marginal_cdf_target(MultiIndex<Real>({1})),
      23           0 :     _conditional_cdf_energy(MultiIndex<Real>({1}))
      24             : {
      25           0 : }
      26             : 
      27           1 : DiscreteFissionPKAPDF::DiscreteFissionPKAPDF(const std::vector<unsigned int> & ZAID,
      28             :                                              const std::vector<Real> & energies,
      29           1 :                                              const MultiIndex<Real> & probabilities)
      30             :   : DiscretePKAPDFBase(ZAID, energies),
      31           1 :     _marginal_cdf_target(probabilities),
      32           1 :     _conditional_cdf_energy(probabilities)
      33             : {
      34             :   // check the size of _probabilities
      35           1 :   if (probabilities.dim() != 2)
      36           0 :     mooseError("probabilities MultiIndex object has wrong dimensions.");
      37             : 
      38           1 :   MultiIndex<Real>::size_type shape = probabilities.size();
      39           1 :   if (shape[0] != _nZA || shape[1] != _ng)
      40           0 :     mooseError("Size of probabilities is inconsistent with random variable input.");
      41             : 
      42           1 :   precomputeCDF(probabilities);
      43           1 :   readFissionData(ZAID);
      44           1 :   computeMagnitude(probabilities);
      45           1 : }
      46             : 
      47             : void
      48           1 : DiscreteFissionPKAPDF::precomputeCDF(MultiIndex<Real> probabilities)
      49             : {
      50             :   // Step 1: Compute the marginal distribution w.r.t. incident energy of the neutron
      51           1 :   MultiIndex<Real>::size_type shape(1);
      52           1 :   MultiIndex<Real>::size_type index(2);
      53             : 
      54           1 :   shape[0] = _nZA;
      55           1 :   _marginal_cdf_target = MultiIndex<Real>(shape);
      56           3 :   for (unsigned int jZA = 0; jZA < _nZA; ++jZA)
      57             :   {
      58             :     Real integral_value = 0.0;
      59           6 :     for (unsigned int jE = 0; jE < _ng; ++jE)
      60           4 :       integral_value += probabilities({jZA, jE});
      61           2 :     _marginal_cdf_target({jZA}) = integral_value;
      62             :   }
      63             : 
      64             :   // Step 2: Compute cdf
      65           2 :   for (auto target : _marginal_cdf_target)
      66             :   {
      67           2 :     index = target.first;
      68           2 :     index[0] -= 1;
      69           2 :     if (target.first[0] > 0)
      70           1 :       target.second += _marginal_cdf_target(index);
      71             :   }
      72             : 
      73             :   // Step 3: Renormalize to ensure that cdf[-1] == 1
      74           1 :   index[0] = _nZA - 1;
      75           1 :   Real last_value = _marginal_cdf_target(index);
      76           2 :   for (auto target : _marginal_cdf_target)
      77           2 :     target.second /= last_value;
      78             : 
      79             :   // Step 1: Compute the conditional distribution of the target w.r.t the sampled energy
      80           1 :   _conditional_cdf_energy = probabilities;
      81             : 
      82             :   // Step 2: Compute cdf
      83           4 :   for (auto energy : _conditional_cdf_energy)
      84             :   {
      85           4 :     index = energy.first;
      86           4 :     index[1] -= 1;
      87           4 :     if (energy.first[1] > 0)
      88           2 :       energy.second += _conditional_cdf_energy(index);
      89             :   }
      90             : 
      91             :   // Step 3: Renormalize to ensure that cdf[-1] == 1
      92             :   std::vector<Real> last_element_per_target;
      93           3 :   for (unsigned int j = 0; j < _conditional_cdf_energy.size()[0]; ++j)
      94           2 :     last_element_per_target.push_back(_conditional_cdf_energy({j, _ng - 1}));
      95           4 :   for (auto energy : _conditional_cdf_energy)
      96           4 :     energy.second /= last_element_per_target[energy.first[0]];
      97           1 : }
      98             : 
      99             : void
     100      500000 : DiscreteFissionPKAPDF::drawSample(std::vector<MyTRIM_NS::IonBase> & initial_state) const
     101             : {
     102             :   // resize initial_state
     103      500000 :   initial_state.resize(2);
     104             : 
     105             :   // first sample the target ZAID and the energy group of the incident neutron
     106             :   std::vector<unsigned int> sampled_indices;
     107      500000 :   sampled_indices.push_back(sampleHelper(_marginal_cdf_target));
     108      500000 :   sampled_indices.push_back(sampleHelper(_conditional_cdf_energy, sampled_indices));
     109             : 
     110             :   // first we need to compute Z and m from ZAID
     111      500000 :   auto target_Z = MagpieUtils::getZFromZAID(_zaids[sampled_indices[0]]);
     112      500000 :   auto target_A = MagpieUtils::getAFromZAID(_zaids[sampled_indices[0]]);
     113             : 
     114             :   // the real random variables also need to be resampled uniformly
     115             :   // within bin index[j]
     116             :   auto neutron_energy =
     117      500000 :       (_energies[sampled_indices[1]] - _energies[sampled_indices[1] + 1]) * MooseRandom::rand() +
     118      500000 :       _energies[sampled_indices[1] + 1];
     119             : 
     120             :   // pull the right fission yield table and sample Z, A, and energy
     121      500000 :   auto energy = MagpieUtils::determineNeutronType(neutron_energy);
     122      500000 :   auto key = _zaids[sampled_indices[0]];
     123             : 
     124             :   // pull the map based on the neutron energy
     125      500000 :   auto zaid_map = _fission_zaids[energy];
     126             :   auto cdf_map = _fission_cdf[energy];
     127             : 
     128             :   /// check if key exists
     129      500000 :   if (zaid_map.find(key) == zaid_map.end())
     130           0 :     mooseError("Fission product ZAID's not found for this target isotope, ",
     131             :                key,
     132             :                " at this energy, ",
     133             :                energy);
     134             : 
     135      500000 :   if (cdf_map.find(key) == cdf_map.end())
     136           0 :     mooseError("Sum yield not found for this target isotope, ", key, " at this energy, ", energy);
     137             : 
     138             :   // now we know the key1/key2 pair exists
     139      500000 :   auto zaid_products = zaid_map[key];
     140      500000 :   auto cdf_products = cdf_map[key];
     141             : 
     142             :   // calculate the total kinetic energy of the fission products
     143      500000 :   auto total_ke = determineFragmentsEnergy(target_Z, target_A);
     144             : 
     145             :   // sample the fission yield cdf to determine the first fission product
     146      500000 :   auto ffindex = sampleHelper(cdf_products);
     147      500000 :   auto ffzaid = zaid_products[ffindex];
     148      500000 :   initial_state[0]._Z = MagpieUtils::getZFromZAID(ffzaid);
     149      500000 :   initial_state[0]._m = MagpieUtils::getAFromZAID(ffzaid);
     150             : 
     151             :   // sample neutrons per fission
     152      500000 :   auto neutrons_per_fission = sampleNu(energy, key);
     153             : 
     154             :   // mass balance to find Z and A of second fission product
     155      500000 :   initial_state[1]._Z = target_Z - initial_state[0]._Z;
     156      500000 :   initial_state[1]._m = target_A - initial_state[0]._m - neutrons_per_fission;
     157             : 
     158             :   // calculate the kinetic energy of each fission product
     159      500000 :   initial_state[0]._E = total_ke / (1.0 + (initial_state[0]._m / initial_state[1]._m));
     160      500000 :   initial_state[1]._E = (initial_state[0]._m / initial_state[1]._m) * initial_state[0]._E;
     161     1000000 : }
     162             : 
     163             : void
     164           1 : DiscreteFissionPKAPDF::readFissionData(const std::vector<unsigned int> & zaid_list)
     165             : {
     166           1 :   _fission_zaids.resize(MagpieUtils::NET_MAX);
     167           1 :   _fission_cdf.resize(MagpieUtils::NET_MAX);
     168             : 
     169             :   // iterate over all neutron energy types
     170           5 :   for (unsigned int energy = 0; energy < MagpieUtils::NET_MAX; ++energy)
     171             :   {
     172             :     // Dummy maps
     173             :     std::map<unsigned int, std::vector<unsigned int>> zaid_map;
     174             :     std::map<unsigned int, std::vector<Real>> cdf_map;
     175          12 :     for (auto & zaid : zaid_list)
     176             :     {
     177             :       // check if $ENDF_FP_DIR is set
     178           8 :       auto path = std::getenv("ENDF_FP_DIR");
     179           8 :       if (path == NULL)
     180           0 :         mooseError("Set $ENDF_FP_DIR to the directory holding the sum yield data files.");
     181             : 
     182             :       // check if file exists, if not continue
     183             :       std::string filename =
     184          16 :           path + std::to_string(zaid) + "_" + MagpieUtils::neutronEnergyName(energy) + ".txt";
     185           8 :       if (!MooseUtils::checkFileReadable(filename, false, false))
     186             :       {
     187             :         // most isotopes have fission data for High but not for Fast wich is weird
     188             :         // if fast is not available, then try to load High instead
     189           1 :         if (MagpieUtils::neutronEnergyName(energy) == "Fast")
     190             :         {
     191           1 :           std::string alt_filename = path + std::to_string(zaid) + "_" + "High" + ".txt";
     192           1 :           if (!MooseUtils::checkFileReadable(alt_filename, false, false))
     193             :             continue;
     194             :           filename = alt_filename;
     195             :         }
     196             :         else
     197           0 :           continue;
     198             :       }
     199             : 
     200             :       // read the data
     201           8 :       std::ifstream infile(filename.c_str());
     202             :       std::vector<unsigned int> zaid_target;
     203             :       std::vector<Real> fission_probabilities;
     204             :       int aa;
     205             :       Real bb;
     206             :       Real prob_prev = 0.0;
     207       15628 :       while (infile >> aa >> bb)
     208             :       {
     209        7806 :         zaid_target.push_back(aa);
     210        7806 :         fission_probabilities.push_back(bb + prob_prev);
     211             : 
     212             :         // accumulate CDF as it reads from the file
     213        7806 :         prob_prev += bb;
     214             :       }
     215             : 
     216             :       // renormalize to ensure that cdf[-1] == 1
     217           8 :       unsigned int last_index = fission_probabilities.size() - 1;
     218           8 :       Real last_value = fission_probabilities[last_index];
     219        7814 :       for (auto & probability : fission_probabilities)
     220        7806 :         probability /= last_value;
     221             : 
     222             :       // Step 4: Store std::vectors into maps
     223           8 :       zaid_map[zaid] = zaid_target;
     224           8 :       cdf_map[zaid] = fission_probabilities;
     225           8 :     }
     226             : 
     227             :     // store map into vector for given energy type
     228           4 :     _fission_zaids[energy] = zaid_map;
     229             :     _fission_cdf[energy] = cdf_map;
     230             :   }
     231           1 : }
     232             : 
     233             : Real
     234      500000 : DiscreteFissionPKAPDF::determineFragmentsEnergy(unsigned int Z, unsigned int A) const
     235             : {
     236             :   // MyTRIM expects energies to be in eV, hence the 1.0e6 factor
     237      500000 :   return (0.1178 * (Utility::pow<2>(Z) / std::cbrt(A)) + 5.8) * 1.0e6;
     238             : }
     239             : 
     240             : /**
     241             :  * Given the target isotope and neutron energy,
     242             :  * determine the avg. number of neutrons per fission.
     243             :  */
     244             : unsigned int
     245      500000 : DiscreteFissionPKAPDF::sampleNu(MagpieUtils::NeutronEnergyType energy_type, unsigned int zaid) const
     246             : {
     247             :   Real nu_bar;
     248      500000 :   if (zaid == 922350 && energy_type == MagpieUtils::Thermal)
     249             :     nu_bar = 2.4355;
     250      500000 :   else if (zaid == 922330 && energy_type == MagpieUtils::Epithermal)
     251             :     nu_bar = 2.4968;
     252      500000 :   else if (zaid == 902320 && energy_type == MagpieUtils::Epithermal)
     253             :     nu_bar = 2.456;
     254      500000 :   else if (zaid == 942390 && energy_type == MagpieUtils::Thermal)
     255             :     nu_bar = 2.8836;
     256           0 :   else if (zaid == 942390 && energy_type == MagpieUtils::Epithermal)
     257             :     nu_bar = 2.8836;
     258           0 :   else if (zaid == 942400 && energy_type == MagpieUtils::Fast)
     259             :     nu_bar = 3.086;
     260           0 :   else if (zaid == 942410 && energy_type == MagpieUtils::Thermal)
     261             :     nu_bar = 2.4979;
     262           0 :   else if (zaid == 942420 && energy_type == MagpieUtils::Epithermal)
     263             :     nu_bar = 3.189;
     264           0 :   else if (zaid == 952410 && energy_type == MagpieUtils::Thermal)
     265             :     nu_bar = 3.239;
     266           0 :   else if (zaid == 962420 && energy_type == MagpieUtils::Epithermal)
     267             :     nu_bar = 2.529;
     268           0 :   else if (zaid == 962430 && energy_type == MagpieUtils::Thermal)
     269             :     nu_bar = 3.433;
     270           0 :   else if (zaid == 962430 && energy_type == MagpieUtils::Epithermal)
     271             :     nu_bar = 3.433;
     272             :   else
     273             :     nu_bar = 2.5;
     274             : 
     275      500000 :   if (nu_bar <= std::floor(nu_bar) + MooseRandom::rand())
     276       58301 :     return std::floor(nu_bar);
     277             :   else
     278      441699 :     return std::ceil(nu_bar);
     279             : }
     280             : 
     281             : void
     282           1 : DiscreteFissionPKAPDF::computeMagnitude(MultiIndex<Real> probabilities)
     283             : {
     284           1 :   _magnitude = 0.0;
     285           4 :   for (auto it : probabilities)
     286             :   {
     287           4 :     MultiIndex<Real>::size_type index = it.first;
     288           4 :     Real delE = _energies[index[1]] - _energies[index[1] + 1];
     289           4 :     _magnitude += delE * it.second;
     290           4 :   }
     291           1 : }