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 "NeutronicsSpectrumSamplerBase.h"
      10             : #include "MooseMesh.h"
      11             : 
      12             : #ifdef RATTLESNAKE_ENABLED
      13             : #include "IsotopeUtilities.h"
      14             : #endif
      15             : 
      16             : // C++ includes
      17             : #include <sstream>
      18             : #include <algorithm>
      19             : #include <limits>
      20             : 
      21             : InputParameters
      22           0 : NeutronicsSpectrumSamplerBase::validParams()
      23             : {
      24           0 :   InputParameters params = ElementUserObject::validParams();
      25           0 :   params.addRequiredParam<std::vector<std::string>>("target_isotope_names",
      26             :                                                     "The list of target isotope names e.g. U235.");
      27           0 :   params.addRequiredCoupledVar("number_densities", "Number densities for each isotope.");
      28           0 :   params.addRequiredParam<std::vector<Real>>("energy_group_boundaries",
      29             :                                              "The energy group boundaries in units of eV orderd "
      30             :                                              "from high to low [natural neutronics ordering].");
      31           0 :   params.addRequiredParam<unsigned int>(
      32             :       "L", "The maximum order of Legendre terms for expanding the recoil XS in mu_lab.");
      33           0 :   params.addRequiredParam<std::vector<Point>>(
      34             :       "points", "The points where you want to evaluate the variables");
      35           0 :   params.addClassDescription("Radiation Damage user object base class.\n Computes PDFs from "
      36             :                              "neutronics calculations that are used to sample PKAs in BCMC "
      37             :                              "simulations.");
      38           0 :   return params;
      39           0 : }
      40             : 
      41           0 : NeutronicsSpectrumSamplerBase::NeutronicsSpectrumSamplerBase(const InputParameters & parameters)
      42             :   : ElementUserObject(parameters),
      43           0 :     _target_isotope_names(getParam<std::vector<std::string>>("target_isotope_names")),
      44           0 :     _energy_group_boundaries(getParam<std::vector<Real>>("energy_group_boundaries")),
      45           0 :     _I(_target_isotope_names.size()),
      46           0 :     _G(_energy_group_boundaries.size() - 1),
      47           0 :     _L(getParam<unsigned int>("L")),
      48           0 :     _points(getParam<std::vector<Point>>("points")),
      49           0 :     _npoints(_points.size()),
      50           0 :     _qp_is_cached(false)
      51             : {
      52             :   // check input dimensions
      53           0 :   if (coupledComponents("number_densities") != _I)
      54           0 :     mooseError("ZAID and number_densities must have the same length.");
      55             : 
      56             :   // get the number density variables
      57           0 :   _number_densities.resize(_I);
      58           0 :   for (unsigned int i = 0; i < _I; ++i)
      59           0 :     _number_densities[i] = &coupledValue("number_densities", i);
      60             : 
      61           0 :   _zaids.resize(_I);
      62           0 :   for (unsigned int i = 0; i < _I; ++i)
      63             :   {
      64             : #ifdef RATTLESNAKE_ENABLED
      65             :     unsigned int Z, A;
      66             :     // check if isotope names are valid, NOTE: error handling is delegated to Yakxs::Utilities
      67             :     YAKXS::Utility::getAZ(_target_isotope_names[i], A, Z);
      68             :     // convert from name to ZAID
      69             :     _zaids[i] = YAKXS::Utility::getZaid(_target_isotope_names[i]);
      70             : #else
      71           0 :     _zaids[i] = localStringToZaid(_target_isotope_names[i]);
      72             : #endif
      73             :   }
      74             : 
      75           0 :   _owner.resize(_npoints);
      76           0 :   _qp_cache.resize(_npoints);
      77             : 
      78             :   // check energy group ordering
      79           0 :   if (_energy_group_boundaries.size() < 2)
      80           0 :     mooseError("At least 2 boundaries must be provided for energy_group_boundaries");
      81           0 :   for (unsigned int j = 1; j < _energy_group_boundaries.size(); ++j)
      82           0 :     if (_energy_group_boundaries[j] >= _energy_group_boundaries[j - 1])
      83           0 :       mooseError("energy_group_boundaries must be ordered from high to low");
      84           0 : }
      85             : 
      86             : void
      87           0 : NeutronicsSpectrumSamplerBase::execute()
      88             : {
      89           0 :   if (_local_elem_to_contained_points.find(_current_elem) != _local_elem_to_contained_points.end())
      90             :   {
      91           0 :     for (auto & j : _local_elem_to_contained_points[_current_elem])
      92             :     {
      93           0 :       _current_point = j;
      94             :       // NOTE: PKA is evaluated at a point that is not necessarily a qp but material
      95             :       // props only live on qps. This is a cheap and dirty solution evaluating the
      96             :       // PKA at the closest quadrature point.
      97           0 :       if (!_qp_is_cached)
      98             :       {
      99           0 :         Real min_dsq = (_q_point[0] - _points[_current_point]).norm_sq();
     100             :         unsigned int min_qp = 0;
     101           0 :         for (unsigned int qp = 1; qp < _q_point.size(); qp++)
     102             :         {
     103             :           Real dsq = (_q_point[qp] - _points[_current_point]).norm_sq();
     104           0 :           if (dsq < min_dsq)
     105             :           {
     106             :             min_dsq = dsq;
     107             :             min_qp = qp;
     108             :           }
     109             :         }
     110           0 :         _qp_cache[_current_point] = min_qp;
     111             :       }
     112             : 
     113             :       // Set current qp to cached min_qp
     114           0 :       _qp = _qp_cache[_current_point];
     115             : 
     116             :       // call back before computing radiation damage PDF for caching things that
     117             :       // don't change
     118           0 :       preComputeRadiationDamagePDF();
     119             : 
     120             :       // Evalute the radiation damage PDF at min_qp
     121           0 :       for (unsigned int i = 0; i < _I; ++i)
     122           0 :         for (unsigned int g = 0; g < _G; ++g)
     123           0 :           for (unsigned int p = 0; p < _nphi; ++p)
     124           0 :             for (unsigned int q = 0; q < _nmu; ++q)
     125           0 :               _sample_point_data[_current_point]({i, g, p, q}) =
     126           0 :                   computeRadiationDamagePDF(i, g, p, q);
     127             :     }
     128             :   }
     129           0 : }
     130             : 
     131             : void
     132           0 : NeutronicsSpectrumSamplerBase::preComputeRadiationDamagePDF()
     133             : {
     134           0 : }
     135             : 
     136             : void
     137           0 : NeutronicsSpectrumSamplerBase::meshChanged()
     138             : {
     139             :   // make sure that _local_elem_to_contained_points is empty
     140             :   _local_elem_to_contained_points.clear();
     141             : 
     142             :   // Rebuild point_locator & find the right element for each point
     143           0 :   std::unique_ptr<PointLocatorBase> point_locator = _mesh.getPointLocator();
     144           0 :   for (unsigned int j = 0; j < _npoints; j++)
     145             :   {
     146             :     // make sure element is local
     147           0 :     _owner[j] = 0;
     148           0 :     const Elem * candidate_elem = (*point_locator)(_points[j]);
     149           0 :     if (candidate_elem && candidate_elem->processor_id() == processor_id())
     150             :     {
     151           0 :       _owner[j] = processor_id();
     152           0 :       if (_local_elem_to_contained_points.find(candidate_elem) !=
     153             :           _local_elem_to_contained_points.end())
     154           0 :         _local_elem_to_contained_points[candidate_elem].push_back(j);
     155             :       else
     156           0 :         _local_elem_to_contained_points[candidate_elem] = {j};
     157             :     }
     158             :   }
     159             : 
     160             :   // make sure everyone knows who owns which point
     161           0 :   for (unsigned int j = 0; j < _npoints; j++)
     162           0 :     gatherMax(_owner[j]);
     163             : 
     164           0 :   _qp_is_cached = false;
     165           0 : }
     166             : 
     167             : void
     168           0 : NeutronicsSpectrumSamplerBase::initialSetup()
     169             : {
     170             :   // allocate PDF
     171             :   // NOTE: Needs to be delayed to initialize because _nSH is set in derived class
     172           0 :   for (unsigned int j = 0; j < _npoints; ++j)
     173           0 :     _sample_point_data.push_back(MultiIndex<Real>({_I, _G, _nphi, _nmu}));
     174           0 : }
     175             : 
     176             : void
     177           0 : NeutronicsSpectrumSamplerBase::initialize()
     178             : {
     179           0 :   meshChanged();
     180             : 
     181           0 :   for (unsigned j = 0; j < _npoints; ++j)
     182           0 :     for (auto entry : _sample_point_data[j])
     183           0 :       entry.second = 0;
     184           0 : }
     185             : 
     186             : void
     187           0 : NeutronicsSpectrumSamplerBase::finalize()
     188             : {
     189           0 :   if (_mesh.n_processors() > 1)
     190             :   {
     191           0 :     for (unsigned j = 0; j < _npoints; ++j)
     192             :     {
     193           0 :       std::vector<Real> flat_data(_I * _G * _nmu * _nphi);
     194           0 :       if (_owner[j] == processor_id())
     195           0 :         flat_data = _sample_point_data[j].getRawData();
     196             : 
     197           0 :       _communicator.broadcast(flat_data, _owner[j]);
     198             : 
     199           0 :       if (_owner[j] != processor_id())
     200           0 :         _sample_point_data[j] = MultiIndex<Real>({_I, _G, _nphi, _nmu}, flat_data);
     201             :     }
     202             :   }
     203             : 
     204             :   // set _qp_is_cached flag to true. Actually we only need to do this if
     205             :   // it was false but this way we can save an if statement
     206           0 :   _qp_is_cached = true;
     207           0 : }
     208             : 
     209             : void
     210           0 : NeutronicsSpectrumSamplerBase::threadJoin(const UserObject & y)
     211             : {
     212             :   const NeutronicsSpectrumSamplerBase & uo = static_cast<const NeutronicsSpectrumSamplerBase &>(y);
     213           0 :   for (unsigned j = 0; j < _npoints; ++j)
     214           0 :     for (unsigned int i = 0; i < _I; ++i)
     215           0 :       for (unsigned int g = 0; g < _G; ++g)
     216           0 :         for (unsigned int p = 0; p < _nphi; ++p)
     217           0 :           for (unsigned int q = 0; q < _nmu; ++q)
     218           0 :             _sample_point_data[j]({i, g, p, q}) += uo._sample_point_data[j]({i, g, p, q});
     219           0 : }
     220             : 
     221             : MultiIndex<Real>
     222           0 : NeutronicsSpectrumSamplerBase::getPDF(unsigned int point_id) const
     223             : {
     224           0 :   return _sample_point_data[point_id];
     225             : }
     226             : 
     227             : std::vector<unsigned int>
     228           0 : NeutronicsSpectrumSamplerBase::getZAIDs() const
     229             : {
     230           0 :   return _zaids;
     231             : }
     232             : 
     233             : std::vector<Real>
     234           0 : NeutronicsSpectrumSamplerBase::getEnergies() const
     235             : {
     236           0 :   return _energy_group_boundaries;
     237             : }
     238             : 
     239             : bool
     240           0 : NeutronicsSpectrumSamplerBase::hasIsotope(std::string target_isotope) const
     241             : {
     242           0 :   auto it = std::find(_target_isotope_names.begin(), _target_isotope_names.end(), target_isotope);
     243           0 :   return it != _target_isotope_names.end();
     244             : }
     245             : 
     246             : unsigned int
     247           0 : NeutronicsSpectrumSamplerBase::localStringToZaid(std::string s) const
     248             : {
     249           0 :   if (s == "U235")
     250             :     return 922350;
     251           0 :   else if (s == "U238")
     252             :     return 922380;
     253           0 :   else if (s == "Pu238")
     254             :     return 942380;
     255           0 :   else if (s == "Pu239")
     256             :     return 942390;
     257           0 :   else if (s == "Pu240")
     258             :     return 942400;
     259           0 :   mooseError(
     260             :       "Isotope name ", s, " cannot be converted with the localStringToZaid conversion method.");
     261             :   return 0;
     262             : }