mass number matching

include/userobjects/MyTRIMRasterizer.h

@@ -66,8 +66,8 @@ class MyTRIMRasterizer : public ElementUserObject
 
   struct PKAParameters
   {
-    /// masses charges (m_i, Z_i) of the matrix elements
-    std::vector<std::pair<Real, Real>> _mass_charge_pair;
+    /// masses, charges, and mass number tolerances (m_i, Z_i, mtol_i) of the matrix elements
+    std::vector<std::tuple<Real, Real, Real>> _mass_charge_tuple;
 
     /// how many isotopes to we have for each element and which index contains the
     /// first Z match? (only support up to Z=119!)

src/userobjects/MyTRIMRasterizer.C

@@ -83,6 +83,8 @@ validParams<MyTRIMRasterizer>()
                        "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>(
@@ -160,10 +162,18 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters)
   auto trim_M = getParam<std::vector<Real>>("M");
   auto trim_Z = getParam<std::vector<Real>>("Z");
 
+  std::vector<Real> mtol;
+  if (isParamValid("Mtol"))
+    mtol = getParam<std::vector<Real>>("Mtol");
+  else
+    mtol.assign(_nvars, 0.5);
+
   if (trim_M.size() != _nvars)
     mooseError("Parameter 'M' must have as many components as coupled variables.");
   if (trim_Z.size() != _nvars)
     mooseError("Parameter 'Z' must have as many components as coupled variables.");
+  if (mtol.size() != _nvars)
+    mooseError("Parameter 'mtol' must have as many components as coupled variables.");
 
   // error check masses and charges
   for (unsigned int i = 0; i < _nvars; ++i)
@@ -252,15 +262,15 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters)
   // setup invariant PKA generation parameters
   for (auto & nZ : _pka_parameters._index_Z)
     nZ = std::make_pair(0, 0);
-  _pka_parameters._mass_charge_pair.resize(_nvars);
+  _pka_parameters._mass_charge_tuple.resize(_nvars);
   _pka_parameters._recoil_rate_scaling = _trim_parameters.recoil_rate_scaling;
   for (unsigned int i = 0; i < _nvars; ++i)
   {
     const auto Z = _trim_parameters.element_prototypes[i]._Z;
 
     // insert (mass, charge) pair
-    _pka_parameters._mass_charge_pair[i] =
-        std::make_pair(_trim_parameters.element_prototypes[i]._m, Z);
+    _pka_parameters._mass_charge_tuple[i] =
+        std::make_tuple(_trim_parameters.element_prototypes[i]._m, Z, mtol[i]);
 
     // increase the count of elements with the same Z
     auto & index_Z = _pka_parameters._index_Z[Z];

src/userobjects/PKAGeneratorBase.C

@@ -40,17 +40,25 @@ PKAGeneratorBase::ionTag(const MyTRIMRasterizer::PKAParameters & pka_parameters,
   // this function relies on the exact representation of whole numbers in IEEE floating point
   // numbers up to a reasonable upper limit [Z < m < 300]
 
+  const auto & mZ = pka_parameters._mass_charge_tuple;
+
   // element not found in rasterizer table
   if (pka_parameters._index_Z[Z].first == 0)
     return -1;
 
+  // only one isotope of this element is present
   if (pka_parameters._index_Z[Z].first == 1)
-    return pka_parameters._index_Z[Z].second;
+  {
+    auto & t = mZ[pka_parameters._index_Z[Z].second];
+    if (std::abs(m - std::get<0>(t)) < std::get<2>(t))
+      return pka_parameters._index_Z[Z].second;
+    else
+      return -1;
+  }
 
-  const auto & mZ = pka_parameters._mass_charge_pair;
-  // start the search at the firsat matching Z
+  // start the search at the first matching Z
   for (auto i = pka_parameters._index_Z[Z].second; i < mZ.size(); ++i)
-    if (mZ[i].second == Z && mZ[i].first == m)
+    if (std::get<1>(mZ[i]) == Z && std::abs(m - std::get<0>(mZ[i])) < std::get<2>(mZ[i]))
       return i;
 
   // no matching mass (isotope) found for the given Z

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
+  [../]
 []