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             : #include "ElasticRecoil.h"
      11             : #include "MagpieUtils.h"
      12             : #include "Function.h"
      13             : 
      14             : // gsl includes
      15             : #include "gsl/gsl_sf_legendre.h"
      16             : #include "gsl/gsl_integration.h"
      17             : 
      18             : registerMooseObject("MagpieApp", ElasticRecoil);
      19             : 
      20             : InputParameters
      21          72 : ElasticRecoil::validParams()
      22             : {
      23          72 :   InputParameters params = ScatteringRecoilCrossSection::validParams();
      24          72 :   params.addClassDescription(
      25             :       "Computes recoil cross sections for elastic scattering events. Allows output to csv.");
      26         144 :   params.addRequiredParam<FunctionName>("scattering_law",
      27             :                                         "Function representing the scattering law for neutrons");
      28          72 :   return params;
      29           0 : }
      30             : 
      31          36 : ElasticRecoil::ElasticRecoil(const InputParameters & parameters)
      32          36 :   : ScatteringRecoilCrossSection(parameters), _scattering_law(getFunction("scattering_law"))
      33             : {
      34          36 :   if (_scattering_cross_section.size() != 1)
      35           0 :     mooseError("ElasticRecoil only allows a single input for scattering_xs parameter. Multiple "
      36             :                "entries are reserved for inelastic scattering modes.");
      37          36 : }
      38             : 
      39             : Real
      40     9825600 : ElasticRecoil::getLabCosine(Real E, Real T, Real /*Q*/) const
      41             : {
      42     9825600 :   Real mu_C = getCMCosine(E, T);
      43     9825600 :   return std::sqrt((1 - mu_C) / 2);
      44             : }
      45             : 
      46             : Real
      47    19552000 : ElasticRecoil::getCMCosine(Real E, Real T, Real /*Q*/) const
      48             : {
      49             :   mooseAssert(E >= 0, "Negative neutron energy");
      50             :   mooseAssert(T >= 0, "Negative recoil energy");
      51             :   // just need to avoid NaN here, E=0 can occur on input when lower energy
      52             :   // bin boundary is E = 0
      53    19552000 :   if (E == 0)
      54             :     return -1;
      55             :   // ensure proper boundaries for CM cosine
      56    19552000 :   Real mu = 1 - 2 * T / E / _gamma;
      57    19552000 :   if (mu < -1)
      58       52800 :     return -1;
      59             :   return mu;
      60             : }
      61             : 
      62             : void
      63           8 : ElasticRecoil::execute()
      64             : {
      65             :   // Loop over all Legendre orders to calculate the coefficients
      66          40 :   for (unsigned int l = 0; l < _L + 1; ++l)
      67             :   {
      68             :     // Loop over all the neutron energy groups
      69          64 :     for (unsigned int g = 0; g < _G; ++g)
      70             :     {
      71             :       // Gets the upper and lower energy values of that group
      72          32 :       Real E_u = _neutron_energy_limits[g];
      73          32 :       Real E_l = _neutron_energy_limits[g + 1];
      74             : 
      75             :       // Loop over all the neutron energies within group g
      76       12832 :       for (unsigned int i_E = 0; i_E < _quad_points.size(); ++i_E)
      77             :       {
      78             :         // Get the energies in Gaussian quadrature points and their respective weights
      79       12800 :         Real E = 0.5 * (E_u - E_l) * _quad_points[i_E] + 0.5 * (E_u + E_l);
      80       12800 :         Real w_E = 0.5 * _quad_weights[i_E] * (E_u - E_l);
      81             : 
      82             :         // Calculate maximum amount of energy transferred to the recoil atom
      83       12800 :         Real T_max = _gamma * E;
      84             : 
      85             :         // Loop over all the possible recoil energy bins
      86       62400 :         for (unsigned int t = 0; t < _T; ++t)
      87             :         {
      88             :           // Gets the upper and lower energy values of that bin
      89       49600 :           Real T_u = _recoil_energy_limits[t];
      90       49600 :           Real T_l = _recoil_energy_limits[t + 1];
      91             : 
      92             :           /*
      93             :            * Calculate possible range of angles according to neutron energy group
      94             :            * and recoil energy bin. This approach avoids unphysical behavior (negative
      95             :            * values due to convergence issue of expansion) of recoil cross section.
      96             :            * The subscript L corresponds to the Lab frame.
      97             :            */
      98       49600 :           Real mu_L_min = getLabCosine(E_u, T_l);
      99       49600 :           Real mu_L_max = getLabCosine(E_l, T_u);
     100             : 
     101             :           // Save maximum and mininum Lab frame cosine values
     102       49600 :           _mu_L[t][g][0] = mu_L_min;
     103       49600 :           _mu_L[t][g][1] = mu_L_max;
     104             : 
     105             :           /*
     106             :            * Elastic scaterring case III: T_max < T_l, stop
     107             :            * When the maximum recoil energy is lower than the lowest energy of the recoil energy bin
     108             :            * t
     109             :            */
     110       49600 :           if (T_max < T_l)
     111       25284 :             continue;
     112             : 
     113             :           /*
     114             :            * Elastic scaterring case II: T_max inside T bin, refit
     115             :            * When the maximum recoil energy is whitin the recoil energy bin t, we need to refit
     116             :            * the quadrature rule for the new interval between T_l and T_max
     117             :            */
     118       24316 :           if (T_max > T_l && T_max < T_u)
     119             :             T_u = T_max;
     120             : 
     121             :           /*
     122             :            * Elastic scaterring case I: T_max > interval, ok
     123             :            * When the maximum recoil energy is greater than the highest energy of the recoil energy
     124             :            * bin t Case II also utilizes this piece of code with T_u = T_max
     125             :            */
     126     9750716 :           for (unsigned int i_T = 0; i_T < _quad_points.size(); ++i_T)
     127             :           {
     128             :             // Get the energies in Gaussian quadrature points and their respective weights
     129     9726400 :             Real T = 0.5 * (T_u - T_l) * _quad_points[i_T] + 0.5 * (T_u + T_l);
     130             : 
     131     9726400 :             Real w_T = 0.5 * _quad_weights[i_T] * (T_u - T_l);
     132             : 
     133             :             // Calculate cosine of recoil angle in the CM frame given the neutron and recoil atom
     134             :             // energies
     135     9726400 :             Real mu_C = getCMCosine(E, T);
     136             : 
     137             :             // Calculate cosine of recoil angle in the Lab frame according to geometry rules
     138     9726400 :             Real mu_L = getLabCosine(E, T);
     139             : 
     140             :             /*
     141             :              * Calculate contribution to cross section coefficients
     142             :              * mu_L is scaled from its possible range of values [mu_L_min, mu_L_max] to fit the
     143             :              * interval [-1,1] of the Legendre polynomials
     144             :              */
     145     9726400 :             Real scaled_mu_L = 2 * (mu_L - mu_L_min) / (mu_L_max - mu_L_min) - 1;
     146     9726400 :             _recoil_coefficient[l][t][g] +=
     147     9726400 :                 1 / _xi_g[g] * _scattering_cross_section[0]->value(E, Point()) *
     148    19452800 :                 _neutron_spectrum.value(E, Point()) * _scattering_law.value(mu_C, Point()) * 2.0 /
     149     9726400 :                 _gamma / E * gsl_sf_legendre_Pl(l, scaled_mu_L) * w_T * w_E;
     150             :           }
     151             :         }
     152             :       }
     153             :     }
     154             :   }
     155           8 : }
     156             : 
     157             : #endif