mass number matching (idaholab#300)

include/userobjects/PKAGeneratorBase.h

@@ -58,6 +58,9 @@ class PKAGeneratorBase : public DiscreteElementUserObject
              const std::vector<Real> & rasterizer_m,
              Real Z,
              Real m) const;
+
+  /// tolerance for mass number matching during tagging
+  const Real _mass_number_tolerance;
 };
 
 #endif // PKAGENERATORBASE_H

src/userobjects/PKAFissionFragmentEmpirical.C

@@ -78,5 +78,6 @@ PKAFissionFragmentEmpirical::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_li
     // add PKAs to list
     ion_list.push_back(ion1);
     ion_list.push_back(ion2);
+std::cout << ion1._Z << " " << ion1._m << " " << ion1._tag << " : " << ion2._Z << " " << ion2._m << " " << ion2._tag << std::endl;
   }
 }

src/userobjects/PKAGeneratorBase.C

@@ -10,16 +10,22 @@
 #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.addParam<Real>("mass_number_tolerance",
+                        0.5,
+                        "Tolerance for matching mass numbers of recoils with rasterizer isotopes.");
+  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),
+    _mass_number_tolerance(getParam<Real>("mass_number_tolerance"))
 {
   setRandomResetFrequency(EXEC_TIMESTEP_END);
 }
@@ -31,22 +37,16 @@ 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,
+                         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)
+      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
+    #mass_number_tolerance = 4
+    relative_density = 1
+  [../]
+  [./rasterizer]
+    type = MyTRIMRasterizer
+    var = 'cC cSi cXe'
+    M =   '12 28  136'
+    Z =   '6  14  54'
+    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/ffe/mass_number_tolerance=4"
+    recover = false
+  [../]
 []