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             : #ifdef GSL_ENABLED
      10             : 
      11             : #include "FunctionDPAUserObject.h"
      12             : #include "PolyatomicDisplacementFunction.h"
      13             : 
      14             : registerMooseObject("MagpieApp", FunctionDPAUserObject);
      15             : 
      16             : InputParameters
      17          24 : FunctionDPAUserObject::validParams()
      18             : {
      19          24 :   InputParameters params = DPAUserObjectBase::validParams();
      20          48 :   params.addRequiredParam<std::vector<std::vector<FunctionName>>>(
      21             :       "damage_functions", "Damage functions for each combinations of projectiles and targets.");
      22          48 :   params.addParam<Real>(
      23             :       "max_energy_step_size",
      24          48 :       100,
      25             :       "The maximum energy step size used for integration and interpolation. Default is 100 eV.");
      26          24 :   params.addClassDescription("Computes the dose in dpa from composition, cross section, damage "
      27             :                              "type, and neutron flux. The damage functions are provided as MOOSE "
      28             :                              "functions where energy is provided via the time arg slot.");
      29          24 :   return params;
      30           0 : }
      31             : 
      32          12 : FunctionDPAUserObject::FunctionDPAUserObject(const InputParameters & parameters)
      33          24 :   : DPAUserObjectBase(parameters), _max_delta_E(getParam<Real>("max_energy_step_size"))
      34             : {
      35             :   // get the damage functions
      36             :   std::vector<std::vector<FunctionName>> nm =
      37          36 :       getParam<std::vector<std::vector<FunctionName>>>("damage_functions");
      38          12 :   _damage_functions.resize(nm.size());
      39          28 :   for (unsigned int j = 0; j < nm.size(); ++j)
      40             :   {
      41          16 :     _damage_functions[j].resize(nm[j].size());
      42          40 :     for (unsigned int i = 0; i < nm[j].size(); ++i)
      43          24 :       _damage_functions[j][i] = &getFunctionByName(nm[j][i]);
      44             :   }
      45          12 : }
      46             : 
      47             : void
      48          12 : FunctionDPAUserObject::initialSetup()
      49             : {
      50          12 :   prepare();
      51             : 
      52             :   // make sure the dimension of _damage_functions is correct
      53          12 :   if (_damage_functions.size() != _atomic_numbers.size())
      54           0 :     paramError("damage_functions",
      55             :                "Size of leading dimension ",
      56             :                _damage_functions.size(),
      57             :                " is not identical to number of isotopes ",
      58             :                _atomic_numbers.size());
      59          28 :   for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
      60          16 :     if (_damage_functions[j].size() != _atomic_numbers.size())
      61           0 :       paramError("damage_functions",
      62             :                  "Size of entry ",
      63             :                  j,
      64             :                  " is ",
      65             :                  _damage_functions[j].size(),
      66             :                  " which is not identical to number of isotopes ",
      67             :                  _atomic_numbers.size());
      68          12 : }
      69             : 
      70             : void
      71          12 : FunctionDPAUserObject::execute()
      72             : {
      73          12 :   accumulateDamage();
      74          12 : }
      75             : 
      76             : void
      77          12 : FunctionDPAUserObject::onCompositionChanged()
      78             : {
      79             :   // compute the integral damage functions
      80          12 :   _integral_damage_functions.resize(_atomic_numbers.size());
      81          28 :   for (unsigned int i = 0; i < _atomic_numbers.size(); ++i)
      82             :   {
      83             :     // resize integral damage function
      84          16 :     _integral_damage_functions[i].resize(_atomic_numbers.size());
      85             : 
      86          40 :     for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
      87             :     {
      88             :       // loop over all energies
      89          24 :       Real energy = 0;
      90          24 :       Real max_energy = getMaxEnergy();
      91          24 :       std::vector<Real> energies(1);
      92          24 :       std::vector<Real> integral(1);
      93             :       // this is the energy time step, initially it must be
      94             :       // small but then can up to _max_delta_E
      95             :       Real deltaE = 1e-5;
      96             :       bool keep_going = true;
      97      217752 :       while (keep_going)
      98             :       {
      99      217728 :         Real energy_old = energy;
     100      217728 :         energy += deltaE;
     101      217728 :         if (energy >= max_energy)
     102             :         {
     103          24 :           energy = max_energy;
     104             :           keep_going = false;
     105             :         }
     106             : 
     107             :         // incremental integration from energy_old to energy
     108             :         std::vector<Real> points;
     109             :         std::vector<Real> weights;
     110      217728 :         PolyatomicDisplacementFunction::gslQuadRule(100, points, weights);
     111             : 
     112             :         Real tint = 0;
     113    21990528 :         for (unsigned int qp = 0; qp < points.size(); ++qp)
     114             :         {
     115    21772800 :           Real e = 0.5 * (points[qp] + 1) * deltaE + energy_old;
     116    21772800 :           Real w = 0.5 * weights[qp] * deltaE;
     117    21772800 :           tint += w * _damage_functions[i][j]->value(e, Point());
     118             :         }
     119             : 
     120             :         // update energies and integral vectors
     121      217728 :         energies.push_back(energy);
     122      217728 :         Real old_integral = integral.back();
     123      217728 :         integral.push_back(tint + old_integral);
     124             : 
     125             :         // now do an adjustment of the timestep but only up to _max_delta_E
     126             :         // the factor of 1.01 was found to be well converged for verification_2.i
     127      217728 :         Real proposed_deltaE = deltaE * 1.01;
     128      217728 :         if (proposed_deltaE < _max_delta_E)
     129             :           deltaE = proposed_deltaE;
     130             :       }
     131             : 
     132             :       // now the linear interpolation object is constructed and stored in integral damage function
     133             :       // var
     134          24 :       _integral_damage_functions[i][j] = LinearInterpolation(energies, integral);
     135             :     }
     136             :   }
     137          12 : }
     138             : 
     139             : Real
     140       32960 : FunctionDPAUserObject::integralDamageFunction(Real T, unsigned int i, unsigned int j) const
     141             : {
     142       32960 :   return _integral_damage_functions[i][j].sample(T);
     143             : }
     144             : 
     145             : void
     146          12 : FunctionDPAUserObject::finalize()
     147             : {
     148          12 : }
     149             : 
     150             : #endif