diff --git a/include/userobjects/MyTRIMRasterizer.h b/include/userobjects/MyTRIMRasterizer.h
index 88a51ec3..727af017 100644
--- a/include/userobjects/MyTRIMRasterizer.h
+++ b/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!)
diff --git a/src/userobjects/MyTRIMRasterizer.C b/src/userobjects/MyTRIMRasterizer.C
index 0877a813..552cb4d2 100644
--- a/src/userobjects/MyTRIMRasterizer.C
+++ b/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];
diff --git a/src/userobjects/PKAGeneratorBase.C b/src/userobjects/PKAGeneratorBase.C
index bbf99f74..a94ad6d4 100644
--- a/src/userobjects/PKAGeneratorBase.C
+++ b/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
diff --git a/tests/polyatomic_cascade/gold/tagging_all_Xe_out.e b/tests/polyatomic_cascade/gold/tagging_all_Xe_out.e
new file mode 100644
index 00000000..177c0e00
Binary files /dev/null and b/tests/polyatomic_cascade/gold/tagging_all_Xe_out.e differ
diff --git a/tests/polyatomic_cascade/gold/tagging_out.e b/tests/polyatomic_cascade/gold/tagging_out.e
new file mode 100644
index 00000000..67410076
Binary files /dev/null and b/tests/polyatomic_cascade/gold/tagging_out.e differ
diff --git a/tests/polyatomic_cascade/tagging.i b/tests/polyatomic_cascade/tagging.i
new file mode 100644
index 00000000..8d92c155
--- /dev/null
+++ b/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
+[]
diff --git a/tests/polyatomic_cascade/tests b/tests/polyatomic_cascade/tests
index 25222cce..97ea75c9 100644
--- a/tests/polyatomic_cascade/tests
+++ b/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
+ [../]
[]