mass number matching (#300)

include/userobjects/MyTRIMRasterizer.h

@@ -58,13 +58,14 @@ class MyTRIMRasterizer : public ElementUserObject
   struct AveragedData
   {
     AveragedData(unsigned int nvars = 0)
-      : _elements(nvars, 0.0), _Z(nvars, 0.0), _M(nvars, 0.0), _site_volume(0.0)
+      : _elements(nvars, 0.0), _Z(nvars, 0.0), _M(nvars, 0.0), _M_tol(nvars, 0.5), _site_volume(0.0)
     {
     }
 
     std::vector<Real> _elements;
     std::vector<Real> _Z;
     std::vector<Real> _M;
+    std::vector<Real> _M_tol;
     Real _site_volume;
   };
 
@@ -86,6 +87,9 @@ class MyTRIMRasterizer : public ElementUserObject
     // Element prototype data (charge, mass, displacement and binding energies)
     std::vector<MyTRIM_NS::Element> element_prototypes;
 
+    /// tolerance for matching mass numbers for tagging
+    std::vector<Real> mtol;
+
     /// conversion factor from chosen length unit to Angstrom
     Real length_scale;
 

include/userobjects/PKAGeneratorBase.h

@@ -56,6 +56,7 @@ class PKAGeneratorBase : public DiscreteElementUserObject
   /// finds the right ion tag; -1 means that the nuclide is not tracked, otherwise the index in the rasterizer nuclide vector must be retrieved
   int ionTag(const std::vector<Real> & rasterizer_Z,
              const std::vector<Real> & rasterizer_m,
+             const std::vector<Real> & mass_number_tolerance,
              Real Z,
              Real m) const;
 };

src/userobjects/MyTRIMRasterizer.C

@@ -76,6 +76,8 @@ InputParameters validParams<MyTRIMRasterizer>()
   params.addCoupledVar("periodic_var", "Optional variables that determines the periodicity. If not supplied the first argument of 'var' will be used.");
   params.addRequiredParam<std::vector<Real>>("M", "Element mass in amu");
   params.addRequiredParam<std::vector<Real>>("Z", "Nuclear charge in e");
+  params.addParam<std::vector<Real>>("Mtol",
+                                     "Tolerance on mass number for tagging PKAs with var id.");
   params.addParam<std::vector<Real>>("Ebind", "Binding energy in eV");
   params.addParam<std::vector<Real>>("Edisp", "Displacement threshold in eV");
   params.addRequiredParam<MaterialPropertyName>("site_volume", "Lattice site volume in nm^3 (regardless of the chosen mesh units)");
@@ -141,6 +143,16 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) :
   if (trim_Z.size() != _nvars)
     mooseError("Parameter 'Z' must have as many components as coupled variables.");
 
+  if (isParamValid("Mtol"))
+  {
+    _trim_parameters.mtol = getParam<std::vector<Real>>("Mtol");
+    if (_trim_parameters.mtol.size() != _nvars)
+      mooseError(
+          "Parameter 'Mtol' (if provided) must have as many components as coupled variables.");
+  }
+  else
+    _trim_parameters.mtol.assign(_nvars, 0.5);
+
   for (unsigned int i = 0; i < _nvars; ++i)
       if (trim_Z[i] > trim_M[i])
         mooseError("Value of Z is larger than value of M for entry ", i);
@@ -265,6 +277,7 @@ MyTRIMRasterizer::execute()
       average._elements[i] += qpvol * (*_var[i])[qp];
       average._M[i] = _trim_parameters.element_prototypes[i]._m;
       average._Z[i] = _trim_parameters.element_prototypes[i]._Z;
+      average._M_tol[i] = _trim_parameters.mtol[i];
     }
 
     // average site volume property

src/userobjects/PKAEmpiricalBase.C

@@ -30,7 +30,7 @@ PKAEmpiricalBase::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_list, Real dt
   const Real m = getM();
   const Real E = getE();
 
-  int tag = ionTag(averaged_data._Z, averaged_data._M, Z, m);
+  int tag = ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, Z, m);
 
   unsigned int num_pka = std::floor(recoil_rate_scaling * dt * vol * getPKARate() + getRandomReal());
 

src/userobjects/PKAFissionFragmentEmpirical.C

@@ -64,8 +64,8 @@ PKAFissionFragmentEmpirical::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_li
 
     // the tag is the element this PKA get registered as upon stopping
     // -1 means the PKA will be ignored
-    ion1._tag = ionTag(averaged_data._Z, averaged_data._M, ion1._Z, ion1._m);
-    ion2._tag = ionTag(averaged_data._Z, averaged_data._M, ion2._Z, ion2._m);
+    ion1._tag = ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, ion1._Z, ion1._m);
+    ion2._tag = ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, ion2._Z, ion2._m);
 
     // set location of the fission event
     setPosition(ion1);

src/userobjects/PKAFissionFragmentNeutronics.C

@@ -59,8 +59,10 @@ PKAFissionFragmentNeutronics::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_l
 
       // the tag is the element this PKA get registered as upon stopping
       // -1 means the PKA will be ignored
-      ion[0]._tag = ionTag(averaged_data._Z, averaged_data._M, ion[0]._Z, ion[0]._m);
-      ion[1]._tag = ionTag(averaged_data._Z, averaged_data._M, ion[1]._Z, ion[1]._m);
+      ion[0]._tag =
+          ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, ion[0]._Z, ion[0]._m);
+      ion[1]._tag =
+          ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, ion[1]._Z, ion[1]._m);
 
       // set location of the fission event
       setPosition(ion[0]);

src/userobjects/PKAFixedPointGenerator.C

@@ -59,7 +59,8 @@ PKAFixedPointGenerator::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_list, R
 
     // the tag is the element this PKA get registered as upon stopping
     // -1 means the PKA will be ignored
-    pka._tag = ionTag(averaged_data._Z, averaged_data._M, pka._Z, pka._m);;
+    pka._tag = ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, pka._Z, pka._m);
+    ;
 
     // set stopping criteria
     pka.setEf();

src/userobjects/PKAGeneratorAlphaDecay.C

@@ -155,8 +155,10 @@ PKAGeneratorAlphaDecay::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_list, R
         ion[1]._Z = parent_Z - 2;
         ion[1]._m = parent_A - 4;
 
-        ion[0]._tag = ionTag(averaged_data._Z, averaged_data._M, ion[0]._Z, ion[0]._m);
-        ion[1]._tag = ionTag(averaged_data._Z, averaged_data._M, ion[1]._Z, ion[1]._m);
+        ion[0]._tag =
+            ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, ion[0]._Z, ion[0]._m);
+        ion[1]._tag =
+            ionTag(averaged_data._Z, averaged_data._M, averaged_data._M_tol, ion[1]._Z, ion[1]._m);
 
         // set the energies; use momentum conservation
         ion[0]._E = decay[l]._alpha_energies;

src/userobjects/PKAGeneratorBase.C

@@ -10,16 +10,18 @@
 #include "MagpieUtils.h"
 #include <algorithm>
 
-template<>
-InputParameters validParams<PKAGeneratorBase>()
+template <>
+InputParameters
+validParams<PKAGeneratorBase>()
 {
   InputParameters params = validParams<DiscreteElementUserObject>();
-  params.addClassDescription("PKA generator user object base class.\n Takes pdf and samples PKAs due to various interactions.");
+  params.addClassDescription("PKA generator user object base class.\n Takes pdf and samples PKAs "
+                             "due to various interactions.");
   return params;
 }
 
-PKAGeneratorBase::PKAGeneratorBase(const InputParameters & parameters) :
-    DiscreteElementUserObject(parameters)
+PKAGeneratorBase::PKAGeneratorBase(const InputParameters & parameters)
+  : DiscreteElementUserObject(parameters)
 {
   setRandomResetFrequency(EXEC_TIMESTEP_END);
 }
@@ -31,22 +33,17 @@ PKAGeneratorBase::setPosition(MyTRIM_NS::IonBase & ion) const
 }
 
 int
-PKAGeneratorBase::ionTag(const std::vector<Real> & rasterizer_Z, const std::vector<Real> & rasterizer_m, Real Z, Real m) const
+PKAGeneratorBase::ionTag(const std::vector<Real> & rasterizer_Z,
+                         const std::vector<Real> & rasterizer_m,
+                         const std::vector<Real> & mass_number_tolerance,
+                         Real Z,
+                         Real m) const
 {
-  // this function relies on the exact representation of whole numbers in IEEE floating point numbers
-  // up to a reasonable upper limit [Z < m < 300]
-  unsigned int count = std::count(rasterizer_Z.begin(), rasterizer_Z.end(), Z);
-  if (count == 1)
-  {
-    const auto & it = std::find(rasterizer_Z.begin(), rasterizer_Z.end(), Z);
-    mooseAssert(it != rasterizer_Z.end(), "Z position in rasterizer_Z was unexpectedly not found");
-    unsigned int index = std::distance(rasterizer_Z.begin(), it);
-    return index;
-  }
-  else if (count > 1)
-    for (unsigned int j = 0; j < rasterizer_Z.size(); ++j)
-      if (rasterizer_Z[j] == Z && rasterizer_m[j] == m)
-        return j;
+  // this function relies on the exact representation of whole numbers in IEEE floating point
+  // numbers up to a reasonable upper limit [Z < m < 300]
+  for (unsigned int j = 0; j < rasterizer_Z.size(); ++j)
+    if (rasterizer_Z[j] == Z && std::abs(rasterizer_m[j] - m) < mass_number_tolerance[j])
+      return j;
   return -1;
 }
 

tests/polyatomic_cascade/tagging.i

@@ -0,0 +1,103 @@
+[Mesh]
+  type = MyTRIMMesh
+  dim = 3
+  nx = 2
+  ny = 2
+  nz = 2
+  xmin = 0
+  xmax = 100
+  ymin = 0
+  ymax = 100
+  zmin = 0
+  zmax = 100
+[]
+
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [./dummy]
+  [../]
+[]
+
+[AuxVariables]
+  [./cC]
+    order = CONSTANT
+    family = MONOMIAL
+    initial_condition = 0.5
+  [../]
+  [./cSi]
+    order = CONSTANT
+    family = MONOMIAL
+    initial_condition = 0.5
+  [../]
+  [./cXe]
+    order = CONSTANT
+    family = MONOMIAL
+    initial_condition = 0.0
+  [../]
+
+  [./int]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+[]
+
+[AuxKernels]
+  [./int]
+    variable = int
+    type = MyTRIMElementResultAux
+    runner = runner
+    ivar = 2
+    defect = INT
+  [../]
+[]
+
+[BCs]
+  [./Periodic]
+    [./all]
+      auto_direction = 'x y z'
+    [../]
+  [../]
+[]
+
+[UserObjects]
+  [./ffe]
+    type = PKAFissionFragmentEmpirical
+    fission_rate = 0.00001
+    relative_density = 1
+  [../]
+  [./rasterizer]
+    type = MyTRIMRasterizer
+    var = 'cC cSi cXe'
+    M =   '12 28  136'
+    Z =   '6  14  54'
+    # MTol = '0.5 0.5 4'
+    site_volume = 0.0404
+    pka_generator = ffe
+    periodic_var = dummy
+    analytical_energy_cutoff = 1e5
+  [../]
+  [./runner]
+    type = MyTRIMElementRun
+    rasterizer = rasterizer
+  [../]
+[]
+
+[Postprocessors]
+  [./total_vacancies]
+    type = ElementIntegralVariablePostprocessor
+    variable = int
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 1
+  nl_abs_tol = 1e-10
+[]
+
+[Outputs]
+  exodus = true
+[]

tests/polyatomic_cascade/tests

@@ -38,4 +38,19 @@
     recover = false
     cli_args = "Outputs/file_base=polyatomic_cascade_recoil_scaled_out UserObjects/constant/pka_rate=0.01 UserObjects/rasterizer/max_pka_count=1000 UserObjects/rasterizer/recoil_rate_scaling=0.8"
   [../]
+  [./tagging_default]
+    type = 'Exodiff'
+    input = 'tagging.i'
+    exodiff = 'tagging_out.e'
+    abs_zero = 1e-6
+    recover = false
+  [../]
+  [./tagging_all_Xe]
+    type = 'Exodiff'
+    input = 'tagging.i'
+    exodiff = 'tagging_all_Xe_out.e'
+    abs_zero = 1e-6
+    cli_args = "Outputs/file_base=tagging_all_Xe_out UserObjects/rasterizer/Mtol='0.5 0.5 4'"
+    recover = false
+  [../]
 []