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 "InelasticRecoil.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", InelasticRecoil);
      19             : 
      20             : InputParameters
      21          24 : InelasticRecoil::validParams()
      22             : {
      23          24 :   InputParameters params = ScatteringRecoilCrossSection::validParams();
      24          24 :   params.addClassDescription(
      25             :       "Computes recoil cross sections for inelastic scattering events. Allows output to csv.");
      26          48 :   params.addRequiredParam<std::vector<Real>>("Q",
      27             :                                              "The Q values for all inelastic reaction channels.");
      28          24 :   return params;
      29           0 : }
      30             : 
      31          12 : InelasticRecoil::InelasticRecoil(const InputParameters & parameters)
      32             :   : ScatteringRecoilCrossSection(parameters),
      33          24 :     _q_values(getParam<std::vector<Real>>("Q")),
      34          12 :     _n_levels(_q_values.size())
      35             : {
      36          12 :   if (_n_levels != _scattering_cross_section.size())
      37           0 :     mooseError("Number of scattering_xs entries must be equal to number of Q values. Each "
      38             :                "inelastic mode is represented by one Q value.");
      39          12 : }
      40             : 
      41             : Real
      42     5561712 : InelasticRecoil::getLabCosine(Real E, Real T, Real Q) const
      43             : {
      44     5561712 :   Real mu_C = getCMCosine(E, T, Q);
      45             :   /*
      46             :    *  need to avoid NaN here from E < Eth , this can occur on input when
      47             :    * El (lower group bound) is supplied to compute LAB mu_min/max, note that
      48             :    * this routine returns one; this is consistent with the case Q = 0, E = 0 => mu_C = -1
      49             :    *
      50             :    * We also need to guard against the case mu_C > 1 which can occur when mu_L_max is computed
      51             :    * Both P1 = (E_l, T_l) and P2 = (E_l, T_u) can correspond to impermissible states; compare Fig. 3
      52             :    * in writeup. P1 and P2 can be outside of the permissible range of energies given by the green
      53             :    * and blue curves.
      54             :    */
      55     5561712 :   Real Eth = std::abs(Q) * (_atomic_mass + 1.0) / _atomic_mass;
      56     5561712 :   if (E <= Eth || mu_C >= 1)
      57             :     return 1;
      58     5561712 :   Real d = Eth / E;
      59     5561712 :   Real f = std::sqrt(1.0 - mu_C * mu_C) / (1.0 / std::sqrt(1 - d) - mu_C);
      60     5561712 :   return std::sqrt(1.0 / (1.0 + f * f));
      61             : }
      62             : 
      63             : Real
      64     5561712 : InelasticRecoil::getCMCosine(Real E, Real T, Real Q) const
      65             : {
      66             :   mooseAssert(E >= 0, "Negative neutron energy");
      67             :   mooseAssert(T >= 0, "Negative recoil energy");
      68             :   // need to avoid NaN here from E < Eth , this can occur on input when
      69             :   // El (lower group bound) is supplied to compute CM mu_min/max
      70     5561712 :   Real Eth = std::abs(Q) * (_atomic_mass + 1.0) / _atomic_mass;
      71     5561712 :   if (E <= Eth)
      72             :     return -1;
      73     5561712 :   Real d = Eth / E;
      74     5561712 :   Real mu = (1.0 - 2.0 * T / E / _gamma - 0.5 * d) / std::sqrt(1.0 - d);
      75     5561712 :   if (mu < -1)
      76             :     return -1;
      77             :   return mu;
      78             : }
      79             : 
      80             : Real
      81       28800 : InelasticRecoil::getMaxRecoilEnergy(Real E, Real Et)
      82             : {
      83       28800 :   if (E <= Et)
      84             :     return 0;
      85       28800 :   Real delta = Et / E;
      86       28800 :   return 0.5 * _gamma * E * (1.0 - 0.5 * delta + std::sqrt(1 - delta));
      87             : }
      88             : 
      89             : Real
      90       28800 : InelasticRecoil::getMinRecoilEnergy(Real E, Real Et)
      91             : {
      92       28800 :   if (E <= Et)
      93             :     return 0;
      94       28800 :   Real delta = Et / E;
      95       28800 :   return 0.5 * _gamma * E * (1.0 - 0.5 * delta - std::sqrt(1 - delta));
      96             : }
      97             : 
      98             : void
      99           4 : InelasticRecoil::execute()
     100             : {
     101             :   // inelastic scattering is isotropic in CM. To remind us that the scattering scattering law
     102             :   // is present, this variable is created and used
     103             :   const Real scattering_law = 0.5;
     104             : 
     105           8 :   for (unsigned int level = 0; level < _n_levels; ++level)
     106             :   {
     107           4 :     Real Q = _q_values[level];
     108           4 :     Real threshold_energy = std::abs(Q) * (_atomic_mass + 1.0) / _atomic_mass;
     109             : 
     110             :     // Loop over all Legendre orders to calculate the coefficients
     111          16 :     for (unsigned int l = 0; l < _L + 1; ++l)
     112             :     {
     113             :       // Loop over all the neutron energy groups
     114          24 :       for (unsigned int g = 0; g < _G; ++g)
     115             :       {
     116             :         // Gets the upper and lower energy values of that group
     117          12 :         Real E_u = _neutron_energy_limits[g];
     118          12 :         Real E_l = _neutron_energy_limits[g + 1];
     119             : 
     120             :         // we can skip this energy group of maximum energy is below or equal to threshold
     121          12 :         if (E_u <= threshold_energy)
     122           0 :           continue;
     123             : 
     124             :         // Loop over all the neutron energies within group g
     125        4812 :         for (unsigned int i_E = 0; i_E < _quad_points.size(); ++i_E)
     126             :         {
     127             :           // Get the energies in Gaussian quadrature points and their respective weights
     128        4800 :           Real E = 0.5 * (E_u - E_l) * _quad_points[i_E] + 0.5 * (E_u + E_l);
     129        4800 :           Real w_E = 0.5 * _quad_weights[i_E] * (E_u - E_l);
     130             : 
     131             :           // This is the dimensionless neutron energy [comp. writeup]
     132        4800 :           Real delta = threshold_energy / E;
     133             : 
     134             :           // Calculate maximum amount of energy transferred to the recoil atom
     135        4800 :           Real T_min = getMinRecoilEnergy(E, threshold_energy);
     136        4800 :           Real T_max = getMaxRecoilEnergy(E, threshold_energy);
     137             : 
     138             :           // Loop over all the possible recoil energy bins
     139       28800 :           for (unsigned int t = 0; t < _T; ++t)
     140             :           {
     141             :             // Gets the upper and lower energy values of that bin
     142       24000 :             Real T_u = _recoil_energy_limits[t];
     143       24000 :             Real T_l = _recoil_energy_limits[t + 1];
     144             : 
     145             :             // Adjust the upper and lower limit according to inelastic scattering limits
     146             :             // [this is described in the document scattering_recoils.pdf]
     147       24000 :             if (E_l < threshold_energy)
     148             :               E_l = threshold_energy;
     149             : 
     150             :             // both called with E_u is intendend [see scattering_recoils.pdf]
     151       24000 :             T_u = std::min(T_u, getMaxRecoilEnergy(E_u, threshold_energy));
     152       24000 :             T_l = std::max(T_l, getMinRecoilEnergy(E_u, threshold_energy));
     153             : 
     154             :             // if T_u <= T_l this recoil range is impossible
     155       24000 :             if (T_u <= T_l)
     156           0 :               continue;
     157             : 
     158             :             // Compute the range of lab cosines. NOTE: this is more complicated than
     159             :             // for elastic scattering because mu_L assumes a minimum for T = |Q| / (A+1)
     160       24000 :             Real mu_L_max = std::max(getLabCosine(E_l, T_l, Q), getLabCosine(E_l, T_u, Q));
     161             :             Real mu_L_min;
     162       24000 :             Real Tp = std::abs(Q) / (_atomic_mass + 1);
     163       24000 :             if (T_l <= Tp && T_u >= Tp)
     164           0 :               mu_L_min = std::sqrt(threshold_energy / E_u);
     165             :             else
     166             :             {
     167       24000 :               if (T_u < Tp)
     168           0 :                 mu_L_min = getLabCosine(E_u, T_u, Q);
     169       24000 :               else if (T_l > Tp)
     170       24000 :                 mu_L_min = getLabCosine(E_u, T_l, Q);
     171             :               else
     172           0 :                 mooseError("Should never get here");
     173             :             }
     174             : 
     175             :             // Save maximum and mininum lab frame cosine values
     176       24000 :             _mu_L[t][g][0] = mu_L_min;
     177       24000 :             _mu_L[t][g][1] = mu_L_max;
     178             : 
     179     9624000 :             for (unsigned int i_T = 0; i_T < _quad_points.size(); ++i_T)
     180             :             {
     181             :               // Get the energies in Gaussian quadrature points and their respective weights
     182     9600000 :               Real T = 0.5 * (T_u - T_l) * _quad_points[i_T] + 0.5 * (T_u + T_l);
     183             : 
     184             :               // T must still adhere to the limits T_min and T_max!
     185     9600000 :               if (T <= T_min || T >= T_max)
     186     4110288 :                 continue;
     187             : 
     188     5489712 :               Real w_T = 0.5 * _quad_weights[i_T] * (T_u - T_l);
     189             : 
     190             :               // Calculate cosine of recoil angle in the Lab frame according to geometry rules
     191     5489712 :               Real mu_L = getLabCosine(E, T, Q);
     192             : 
     193             :               /*
     194             :                * Calculate contribution to cross section coefficients
     195             :                * mu_L is scaled from its possible range of values [mu_L_min, mu_L_max] to fit the
     196             :                * interval [-1,1] of the Legendre polynomials
     197             :                */
     198     5489712 :               Real scaled_mu_L = 2 * (mu_L - mu_L_min) / (mu_L_max - mu_L_min) - 1;
     199             : 
     200     5489712 :               Real jacobian_determinant = 2.0 / (_gamma * E * std::sqrt(1 - delta));
     201     5489712 :               _recoil_coefficient[l][t][g] +=
     202    10979424 :                   1 / _xi_g[g] * _scattering_cross_section[level]->value(E, Point()) *
     203     5489712 :                   _neutron_spectrum.value(E, Point()) * scattering_law * jacobian_determinant *
     204     5489712 :                   gsl_sf_legendre_Pl(l, scaled_mu_L) * w_T * w_E;
     205             :             }
     206             :           }
     207             :         }
     208             :       }
     209             :     }
     210             :   }
     211           4 : }
     212             : 
     213             : #endif