From b820a1e4ae43b8aa6bc87aace022afd830fe24c4 Mon Sep 17 00:00:00 2001 From: Daniel Schwen Date: Mon, 13 Aug 2018 08:10:51 -0600 Subject: [PATCH 1/3] Remove mass and Z from average data, refactor ionTag() (#333) --- include/userobjects/MyTRIMRasterizer.h | 34 ++++- include/userobjects/PKAEmpiricalBase.h | 8 +- .../userobjects/PKAFissionFragmentEmpirical.h | 8 +- .../PKAFissionFragmentNeutronics.h | 6 +- include/userobjects/PKAFixedPointGenerator.h | 8 +- include/userobjects/PKAGeneratorAlphaDecay.h | 6 +- include/userobjects/PKAGeneratorBase.h | 11 +- include/userobjects/PKAGeneratorRecoil.h | 6 +- src/userobjects/MyTRIMRasterizer.C | 125 ++++++++++++------ src/userobjects/PKAEmpiricalBase.C | 22 ++- src/userobjects/PKAFissionFragmentEmpirical.C | 35 +++-- .../PKAFissionFragmentNeutronics.C | 39 ++++-- src/userobjects/PKAFixedPointGenerator.C | 26 ++-- src/userobjects/PKAGeneratorAlphaDecay.C | 65 +++++---- src/userobjects/PKAGeneratorBase.C | 52 +++++--- src/userobjects/PKAGeneratorRecoil.C | 38 ++++-- src/userobjects/PKAGun.C | 13 +- 17 files changed, 316 insertions(+), 186 deletions(-) diff --git a/include/userobjects/MyTRIMRasterizer.h b/include/userobjects/MyTRIMRasterizer.h index 51eb5fb3..4695be6b 100644 --- a/include/userobjects/MyTRIMRasterizer.h +++ b/include/userobjects/MyTRIMRasterizer.h @@ -12,6 +12,7 @@ #include "ElementUserObject.h" #include +#include #include #include "mytrim/ion.h" @@ -57,17 +58,33 @@ class MyTRIMRasterizer : public ElementUserObject /// element averaged data struct AveragedData { - AveragedData(unsigned int nvars = 0) - : _elements(nvars, 0.0), _Z(nvars, 0.0), _M(nvars, 0.0), _site_volume(0.0) - { - } + AveragedData(unsigned int nvars = 0) : _elements(nvars, 0.0), _site_volume(0.0) {} std::vector _elements; - std::vector _Z; - std::vector _M; Real _site_volume; }; + struct PKAParameters + { + /// masses charges (m_i, Z_i) of the matrix elements + std::vector> _mass_charge_pair; + + /// how many isotopes to we have for each element? (only support up to Z=119!) + std::array _num_Z; + + /// if only one element with a specific Z is present point to its rasterizer index here for fast lookup + std::array _single_Z_index; + + /// time interval over which the PKAs are added + Real _dt; + + /// recoil rate scaling + Real _recoil_rate_scaling; + + /// current element volume + Real _volume; + }; + enum TRIMModuleEnum { MYTRIM_CORE = 0, @@ -141,6 +158,9 @@ class MyTRIMRasterizer : public ElementUserObject /// Simulation parameters TrimParameters _trim_parameters; + /// Global (non-spatial dependent) parameters required for PKA generation + PKAParameters _pka_parameters; + /// coupled variable values std::vector _var; @@ -152,7 +172,7 @@ class MyTRIMRasterizer : public ElementUserObject std::vector _pka_generators; /// @} - /// @{ material map for the TRIM simulation + /// @{ material map for the TRIM simulation (spatial dependent) typedef std::map MaterialMap; MaterialMap _material_map; /// @} diff --git a/include/userobjects/PKAEmpiricalBase.h b/include/userobjects/PKAEmpiricalBase.h index 3f5dca8d..a54c3356 100644 --- a/include/userobjects/PKAEmpiricalBase.h +++ b/include/userobjects/PKAEmpiricalBase.h @@ -24,11 +24,9 @@ class PKAEmpiricalBase : public PKAGeneratorBase public: PKAEmpiricalBase(const InputParameters & parameters); - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, - const MyTRIMRasterizer::AveragedData & averaged_data) const; + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, + const MyTRIMRasterizer::AveragedData &) const; protected: /// Fission rate (per unit volume) diff --git a/include/userobjects/PKAFissionFragmentEmpirical.h b/include/userobjects/PKAFissionFragmentEmpirical.h index f935e8d0..0954fa4d 100644 --- a/include/userobjects/PKAFissionFragmentEmpirical.h +++ b/include/userobjects/PKAFissionFragmentEmpirical.h @@ -25,11 +25,9 @@ class PKAFissionFragmentEmpirical : public PKAGeneratorBase public: PKAFissionFragmentEmpirical(const InputParameters & parameters); - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, - const MyTRIMRasterizer::AveragedData & averaged_data) const; + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, + const MyTRIMRasterizer::AveragedData &) const; protected: /// Fission rate (per unit volume) assuming pure fully dense UO2 diff --git a/include/userobjects/PKAFissionFragmentNeutronics.h b/include/userobjects/PKAFissionFragmentNeutronics.h index c45b7cef..484ac69f 100644 --- a/include/userobjects/PKAFissionFragmentNeutronics.h +++ b/include/userobjects/PKAFissionFragmentNeutronics.h @@ -23,10 +23,8 @@ class PKAFissionFragmentNeutronics : public PKAGeneratorNeutronicsBase public: PKAFissionFragmentNeutronics(const InputParameters & parameters); - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, const MyTRIMRasterizer::AveragedData &) const; virtual void setPDF(const std::vector & ZAID, diff --git a/include/userobjects/PKAFixedPointGenerator.h b/include/userobjects/PKAFixedPointGenerator.h index 89cd7760..ef81e161 100644 --- a/include/userobjects/PKAFixedPointGenerator.h +++ b/include/userobjects/PKAFixedPointGenerator.h @@ -25,11 +25,9 @@ class PKAFixedPointGenerator : public PKAGeneratorBase public: PKAFixedPointGenerator(const InputParameters & parameters); - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, - const MyTRIMRasterizer::AveragedData & averaged_data) const; + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, + const MyTRIMRasterizer::AveragedData &) const; virtual void meshChanged() { updateCachedElementID(); } protected: diff --git a/include/userobjects/PKAGeneratorAlphaDecay.h b/include/userobjects/PKAGeneratorAlphaDecay.h index d8671c78..37c1209a 100644 --- a/include/userobjects/PKAGeneratorAlphaDecay.h +++ b/include/userobjects/PKAGeneratorAlphaDecay.h @@ -24,10 +24,8 @@ class PKAGeneratorAlphaDecay : public PKAGeneratorBase public: PKAGeneratorAlphaDecay(const InputParameters & parameters); - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, const MyTRIMRasterizer::AveragedData &) const override; /// this stores the decay information for a single alpha decaying nuclide diff --git a/include/userobjects/PKAGeneratorBase.h b/include/userobjects/PKAGeneratorBase.h index f1e156a0..89174d5b 100644 --- a/include/userobjects/PKAGeneratorBase.h +++ b/include/userobjects/PKAGeneratorBase.h @@ -31,10 +31,8 @@ class PKAGeneratorBase : public DiscreteElementUserObject * Append the ions for the current element and time window dt. * The element volume is passed in as it is computed in the MyTRIMRasterizer anyways. */ - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, const MyTRIMRasterizer::AveragedData &) const = 0; virtual void initialize() {} @@ -54,10 +52,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 & rasterizer_Z, - const std::vector & rasterizer_m, - Real Z, - Real m) const; + int ionTag(const MyTRIMRasterizer::PKAParameters &, Real Z, Real m) const; }; #endif // PKAGENERATORBASE_H diff --git a/include/userobjects/PKAGeneratorRecoil.h b/include/userobjects/PKAGeneratorRecoil.h index 3ba504ee..ebc475c2 100644 --- a/include/userobjects/PKAGeneratorRecoil.h +++ b/include/userobjects/PKAGeneratorRecoil.h @@ -23,10 +23,8 @@ class PKAGeneratorRecoil : public PKAGeneratorNeutronicsBase public: PKAGeneratorRecoil(const InputParameters & parameters); - virtual void appendPKAs(std::vector & ion_list, - Real dt, - Real vol, - Real recoil_rate_scaling, + virtual void appendPKAs(std::vector &, + const MyTRIMRasterizer::PKAParameters &, const MyTRIMRasterizer::AveragedData &) const; virtual void setPDF(const std::vector & ZAID, diff --git a/src/userobjects/MyTRIMRasterizer.C b/src/userobjects/MyTRIMRasterizer.C index 81fb1ca8..f831735a 100644 --- a/src/userobjects/MyTRIMRasterizer.C +++ b/src/userobjects/MyTRIMRasterizer.C @@ -16,14 +16,14 @@ #include "libmesh/parallel_algebra.h" // custom data load and data store methods for a struct with an std::vector member -template<> +template <> inline void dataStore(std::ostream & stream, MyTRIMRasterizer::AveragedData & ad, void * context) { dataStore(stream, ad._elements, context); dataStore(stream, ad._site_volume, context); } -template<> +template <> inline void dataLoad(std::istream & stream, MyTRIMRasterizer::AveragedData & ad, void * context) { @@ -31,8 +31,9 @@ dataLoad(std::istream & stream, MyTRIMRasterizer::AveragedData & ad, void * cont dataLoad(stream, ad._site_volume, context); } -// custom data load and data store methods for a class with virtual members (vtable pointer must not be (un)serialized) -template<> +// custom data load and data store methods for a class with virtual members (vtable pointer must not +// be (un)serialized) +template <> inline void dataStore(std::ostream & stream, MyTRIM_NS::IonBase & pl, void * context) { @@ -48,7 +49,7 @@ dataStore(std::ostream & stream, MyTRIM_NS::IonBase & pl, void * context) dataStore(stream, pl._Ef, context); dataStore(stream, pl._state, context); } -template<> +template <> inline void dataLoad(std::istream & stream, MyTRIM_NS::IonBase & pl, void * context) { @@ -67,19 +68,25 @@ dataLoad(std::istream & stream, MyTRIM_NS::IonBase & pl, void * context) registerMooseObject("MagpieApp", MyTRIMRasterizer); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addClassDescription("Gather the element distribution of the simulation domain for a TRIM binary collision Monte Carlo simulation"); + params.addClassDescription("Gather the element distribution of the simulation domain for a TRIM " + "binary collision Monte Carlo simulation"); params.addRequiredCoupledVar("var", "Variables to rasterize"); - params.addCoupledVar("periodic_var", "Optional variables that determines the periodicity. If not supplied the first argument of 'var' will be used."); + params.addCoupledVar("periodic_var", + "Optional variables that determines the periodicity. If not supplied the " + "first argument of 'var' will be used."); params.addRequiredParam>("M", "Element mass in amu"); params.addRequiredParam>("Z", "Nuclear charge in e"); params.addParam>("Ebind", "Binding energy in eV"); params.addParam>("Edisp", "Displacement threshold in eV"); - params.addRequiredParam("site_volume", "Lattice site volume in nm^3 (regardless of the chosen mesh units)"); - params.addRequiredParam>("pka_generator", "List of PKA generating user objects"); + params.addRequiredParam( + "site_volume", "Lattice site volume in nm^3 (regardless of the chosen mesh units)"); + params.addRequiredParam>("pka_generator", + "List of PKA generating user objects"); ExecFlagEnum setup_options(MooseUtils::getDefaultExecFlagEnum()); // we run this object once a timestep @@ -88,27 +95,42 @@ InputParameters validParams() // which TRIM Module to run for optional capabilities like energy deposition MooseEnum trim_module_options("CORE=0 ENERGY_DEPOSITION=1", "CORE"); - params.addParam("trim_module", trim_module_options, "TRIM Module to run for optional capabilities like energy deposition"); + params.addParam("trim_module", + trim_module_options, + "TRIM Module to run for optional capabilities like energy deposition"); // which units of length to use MooseEnum length_unit_options("ANGSTROM=0 NANOMETER MICROMETER", "ANGSTROM"); - params.addParam("length_unit", length_unit_options, "Length units of the MOOSE mesh. MyTRIM contains pretabulated crossection data with units so this option must be set correctly to obtain physical results."); + params.addParam("length_unit", + length_unit_options, + "Length units of the MOOSE mesh. MyTRIM contains pretabulated " + "crossection data with units so this option must be set correctly to " + "obtain physical results."); // Advanced options params.addParam("interval", 1, "The time step interval at which TRIM BCMC is run"); - params.addParam("analytical_energy_cutoff", 0.0, "Energy cutoff in eV below which recoils are not followed explicitly but effects are calculated analytically."); - params.addParam("recoil_rate_scaling", 1.0, "A factor to scale computed reaction rates in the the PKAGenerator objects. This is useful to avoid extremely large PKA lists."); - params.addParam("max_pka_count", "Desired number of PKAs to be run during each invocation of mytrim"); + params.addParam("analytical_energy_cutoff", + 0.0, + "Energy cutoff in eV below which recoils are not followed explicitly but " + "effects are calculated analytically."); + params.addParam("recoil_rate_scaling", + 1.0, + "A factor to scale computed reaction rates in the the PKAGenerator " + "objects. This is useful to avoid extremely large PKA lists."); + params.addParam( + "max_pka_count", "Desired number of PKAs to be run during each invocation of mytrim"); params.addParamNamesToGroup("interval analytical_energy_cutoff max_pka_count", "Advanced"); - params.addParam("r_rec", "Recombination radius in Angstrom. Frenkel pairs with a separation distance lower than this will be removed from the cascade"); + params.addParam("r_rec", + "Recombination radius in Angstrom. Frenkel pairs with a separation " + "distance lower than this will be removed from the cascade"); params.addParamNamesToGroup("r_rec", "Recombination"); return params; } -MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) : - ElementUserObject(parameters), +MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) + : ElementUserObject(parameters), _nvars(coupledComponents("var")), _dim(_mesh.dimension()), _var(_nvars), @@ -142,8 +164,8 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) : mooseError("Parameter 'Z' must have as many components as coupled variables."); 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); + if (trim_Z[i] > trim_M[i]) + mooseError("Value of Z is larger than value of M for entry ", i); for (unsigned int i = 0; i < _nvars; ++i) { @@ -160,7 +182,8 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) : for (unsigned int i = 0; i < _nvars; ++i) _trim_parameters.element_prototypes[i]._Elbind = 3.0; else - mooseError("Parameter 'Ebind' must have as many components as coupled variables (or left empty for a default of 3eV)."); + mooseError("Parameter 'Ebind' must have as many components as coupled variables (or left empty " + "for a default of 3eV)."); auto trim_Edisp = getParam>("Edisp"); if (trim_Edisp.size() == _nvars) @@ -170,7 +193,8 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) : for (unsigned int i = 0; i < _nvars; ++i) _trim_parameters.element_prototypes[i]._Edisp = 25.0; else - mooseError("Parameter 'Edisp' must have as many components as coupled variables (or left empty for a default of 25eV)."); + mooseError("Parameter 'Edisp' must have as many components as coupled variables (or left empty " + "for a default of 25eV)."); if (isParamValid("r_rec")) { @@ -214,6 +238,28 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) : default: mooseError("Unknown length unit."); } + + // setup invariant PKA generation parameters + for (auto & nZ : _pka_parameters._num_Z) + nZ = 0; + _pka_parameters._mass_charge_pair.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); + + // increase the count of elements with the same Z + auto & num_Z = _pka_parameters._num_Z[Z]; + num_Z++; + + // only set this to the first index (important for ionTag()) + if (num_Z == 1) + _pka_parameters._single_Z_index[Z] = i; + } } bool @@ -238,6 +284,9 @@ MyTRIMRasterizer::initialize() _material_map.clear(); _pka_list.clear(); } + + /// setup global PKA parameters for the current timestep + _pka_parameters._dt = _accumulated_time + _fe_problem.dt(); } void @@ -260,12 +309,7 @@ MyTRIMRasterizer::execute() // average compositions on the element for (unsigned int i = 0; i < _nvars; ++i) - { - 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 site volume property average._site_volume += qpvol * _site_volume_prop[qp]; @@ -283,13 +327,16 @@ MyTRIMRasterizer::execute() // store in map _material_map[_current_elem->id()] = average; + // update corrent element volume + _pka_parameters._volume = vol; + // add PKAs for current element for (auto && gen : _pka_generators) - gen->appendPKAs(_pka_list, _accumulated_time + _fe_problem.dt(), vol, _trim_parameters.recoil_rate_scaling, average); + gen->appendPKAs(_pka_list, _pka_parameters, average); } void -MyTRIMRasterizer::threadJoin(const UserObject &y) +MyTRIMRasterizer::threadJoin(const UserObject & y) { // if the map needs to be updated we merge the maps from all threads if (_execute_this_timestep) @@ -348,8 +395,7 @@ MyTRIMRasterizer::finalize() // sort PKA list only on processor 0 & assign random number seeds std::sort(_pka_list.begin(), _pka_list.end(), [](MyTRIM_NS::IonBase a, MyTRIM_NS::IonBase b) { - return (a._pos < b._pos) || - (a._pos == b._pos && a._m < b._m) || + return (a._pos < b._pos) || (a._pos == b._pos && a._m < b._m) || (a._pos == b._pos && a._m == b._m && a._E < b._E) || (a._pos == b._pos && a._m == b._m && a._E == b._E && a._Z < b._Z); }); @@ -361,7 +407,8 @@ MyTRIMRasterizer::finalize() // rejection is performed in processor 0 only _trim_parameters.original_npka = _pka_list.size(); - if (_trim_parameters.desired_npka == 0 || _trim_parameters.desired_npka > _trim_parameters.original_npka) + if (_trim_parameters.desired_npka == 0 || + _trim_parameters.desired_npka > _trim_parameters.original_npka) { _trim_parameters.scaled_npka = _trim_parameters.original_npka; _trim_parameters.result_scaling_factor = 1.0 / _trim_parameters.recoil_rate_scaling; @@ -370,7 +417,8 @@ MyTRIMRasterizer::finalize() { if (processor_id() == 0) { - Real acceptance_probability = Real(_trim_parameters.desired_npka) / Real(_trim_parameters.original_npka); + Real acceptance_probability = + Real(_trim_parameters.desired_npka) / Real(_trim_parameters.original_npka); // most straight-forward but probably inefficient implementation of rejection std::vector old_pka_list = _pka_list; @@ -381,7 +429,9 @@ MyTRIMRasterizer::finalize() // save the size of the PKA list after rejection & the scaling factor _trim_parameters.scaled_npka = _pka_list.size(); - _trim_parameters.result_scaling_factor = Real(_trim_parameters.original_npka) / Real(_trim_parameters.scaled_npka) / _trim_parameters.recoil_rate_scaling; + _trim_parameters.result_scaling_factor = Real(_trim_parameters.original_npka) / + Real(_trim_parameters.scaled_npka) / + _trim_parameters.recoil_rate_scaling; } // need to broadcast the size of the PKA list after rejection & result scaling factor @@ -417,7 +467,7 @@ MyTRIMRasterizer::finalize() // split PKAs into per-processor ranges std::vector interval(_app.n_processors() + 1, 0); for (unsigned int i = 0; i < _app.n_processors(); ++i) - interval[i+1] = (_pka_list.size() - interval[i]) / (_app.n_processors() - i) + interval[i]; + interval[i + 1] = (_pka_list.size() - interval[i]) / (_app.n_processors() - i) + interval[i]; auto begin = interval[processor_id()]; auto end = interval[processor_id() + 1]; @@ -486,7 +536,8 @@ MyTRIMRasterizer::serialize(std::string & serialized_buffer) void MyTRIMRasterizer::deserialize(std::vector & serialized_buffers) { - mooseAssert(serialized_buffers.size() == _app.n_processors(), "Unexpected size of serialized_buffers: " << serialized_buffers.size()); + mooseAssert(serialized_buffers.size() == _app.n_processors(), + "Unexpected size of serialized_buffers: " << serialized_buffers.size()); // The input string stream used for deserialization std::istringstream iss; diff --git a/src/userobjects/PKAEmpiricalBase.C b/src/userobjects/PKAEmpiricalBase.C index 455c8d0a..d3e5ac25 100644 --- a/src/userobjects/PKAEmpiricalBase.C +++ b/src/userobjects/PKAEmpiricalBase.C @@ -8,21 +8,28 @@ #include "PKAEmpiricalBase.h" -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); return params; } -PKAEmpiricalBase::PKAEmpiricalBase(const InputParameters & parameters) : - PKAGeneratorBase(parameters) +PKAEmpiricalBase::PKAEmpiricalBase(const InputParameters & parameters) + : PKAGeneratorBase(parameters) { } void -PKAEmpiricalBase::appendPKAs(std::vector & ion_list, Real dt, Real vol, Real recoil_rate_scaling, const MyTRIMRasterizer::AveragedData & averaged_data) const +PKAEmpiricalBase::appendPKAs(std::vector & ion_list, + const MyTRIMRasterizer::PKAParameters & pka_parameters, + const MyTRIMRasterizer::AveragedData &) const { + const auto dt = pka_parameters._dt; + const auto vol = pka_parameters._volume; + const auto recoil_rate_scaling = pka_parameters._recoil_rate_scaling; + mooseAssert(dt >= 0, "Passed a negative time window into PKAEmpiricalBase::appendPKAs"); mooseAssert(vol >= 0, "Passed a negative volume into PKAEmpiricalBase::appendPKAs"); @@ -30,9 +37,10 @@ PKAEmpiricalBase::appendPKAs(std::vector & 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(pka_parameters, Z, m); - unsigned int num_pka = std::floor(recoil_rate_scaling * dt * vol * getPKARate() + getRandomReal()); + unsigned int num_pka = + std::floor(recoil_rate_scaling * dt * vol * getPKARate() + getRandomReal()); for (unsigned i = 0; i < num_pka; ++i) { diff --git a/src/userobjects/PKAFissionFragmentEmpirical.C b/src/userobjects/PKAFissionFragmentEmpirical.C index 456ac784..b6367c9c 100644 --- a/src/userobjects/PKAFissionFragmentEmpirical.C +++ b/src/userobjects/PKAFissionFragmentEmpirical.C @@ -10,29 +10,42 @@ registerMooseObject("MagpieApp", PKAFissionFragmentEmpirical); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addParam("fission_rate", 1e-8, "Fission rate per unit volume (uses mesh units defined in the rasterizer and moose time units)"); - params.addRequiredCoupledVar("relative_density", "Relative UO2 density (1 is fully dense, 0 is no UO2"); + params.addParam("fission_rate", + 1e-8, + "Fission rate per unit volume (uses mesh units defined in the " + "rasterizer and moose time units)"); + params.addRequiredCoupledVar("relative_density", + "Relative UO2 density (1 is fully dense, 0 is no UO2"); return params; } -PKAFissionFragmentEmpirical::PKAFissionFragmentEmpirical(const InputParameters & parameters) : - PKAGeneratorBase(parameters), +PKAFissionFragmentEmpirical::PKAFissionFragmentEmpirical(const InputParameters & parameters) + : PKAGeneratorBase(parameters), _fission_rate(getPostprocessorValue("fission_rate")), _relative_density(coupledValue("relative_density")) { } void -PKAFissionFragmentEmpirical::appendPKAs(std::vector & ion_list, Real dt, Real vol, Real recoil_rate_scaling, const MyTRIMRasterizer::AveragedData & averaged_data) const +PKAFissionFragmentEmpirical::appendPKAs(std::vector & ion_list, + const MyTRIMRasterizer::PKAParameters & pka_parameters, + const MyTRIMRasterizer::AveragedData &) const { - mooseAssert(dt >= 0, "Passed a negative time window into PKAFissionFragmentEmpirical::appendPKAs"); + const auto dt = pka_parameters._dt; + const auto vol = pka_parameters._volume; + const auto recoil_rate_scaling = pka_parameters._recoil_rate_scaling; + + mooseAssert(dt >= 0, + "Passed a negative time window into PKAFissionFragmentEmpirical::appendPKAs"); mooseAssert(vol >= 0, "Passed a negative volume into PKAFissionFragmentEmpirical::appendPKAs"); - unsigned int num_fission = std::floor(recoil_rate_scaling * dt * vol * _fission_rate * _relative_density[0] + getRandomReal()); + unsigned int num_fission = std::floor( + recoil_rate_scaling * dt * vol * _fission_rate * _relative_density[0] + getRandomReal()); /// Mass inverter to sample PKA mass distribution MyTRIM_NS::MassInverter mass_inverter; @@ -64,8 +77,8 @@ PKAFissionFragmentEmpirical::appendPKAs(std::vector & 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(pka_parameters, ion1._Z, ion1._m); + ion2._tag = ionTag(pka_parameters, ion2._Z, ion2._m); // set location of the fission event setPosition(ion1); diff --git a/src/userobjects/PKAFissionFragmentNeutronics.C b/src/userobjects/PKAFissionFragmentNeutronics.C index c6a1a69d..164d0a47 100644 --- a/src/userobjects/PKAFissionFragmentNeutronics.C +++ b/src/userobjects/PKAFissionFragmentNeutronics.C @@ -13,29 +13,40 @@ registerMooseObject("MagpieApp", PKAFissionFragmentNeutronics); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addClassDescription("PKA generator (fission) user object.\n Takes pdf and samples PKAs due to fission."); + params.addClassDescription( + "PKA generator (fission) user object.\n Takes pdf and samples PKAs due to fission."); return params; } -PKAFissionFragmentNeutronics::PKAFissionFragmentNeutronics(const InputParameters & parameters): - PKAGeneratorNeutronicsBase(parameters) +PKAFissionFragmentNeutronics::PKAFissionFragmentNeutronics(const InputParameters & parameters) + : PKAGeneratorNeutronicsBase(parameters) { } void -PKAFissionFragmentNeutronics::setPDF(const std::vector & ZAID, const std::vector & energies, const MultiIndex & probabilities) +PKAFissionFragmentNeutronics::setPDF(const std::vector & ZAID, + const std::vector & energies, + const MultiIndex & probabilities) { _pdf = DiscreteFissionPKAPDF(ZAID, energies, probabilities); } void -PKAFissionFragmentNeutronics::appendPKAs(std::vector & ion_list, Real dt, Real vol, Real recoil_rate_scaling, const MyTRIMRasterizer::AveragedData & averaged_data) const +PKAFissionFragmentNeutronics::appendPKAs(std::vector & ion_list, + const MyTRIMRasterizer::PKAParameters & pka_parameters, + const MyTRIMRasterizer::AveragedData & averaged_data) const { - mooseAssert(dt >= 0, "Passed a negative time window into PKAFissionFragmentNeutronics::appendPKAs"); + const auto dt = pka_parameters._dt; + const auto vol = pka_parameters._volume; + const auto recoil_rate_scaling = pka_parameters._recoil_rate_scaling; + + mooseAssert(dt >= 0, + "Passed a negative time window into PKAFissionFragmentNeutronics::appendPKAs"); mooseAssert(vol >= 0, "Passed a negative volume into PKAFissionFragmentNeutronics::appendPKAs"); if (averaged_data._elements.size() != _partial_neutronics_reaction_rates.size()) @@ -43,9 +54,11 @@ PKAFissionFragmentNeutronics::appendPKAs(std::vector & ion_l for (unsigned int nuclide = 0; nuclide < _partial_neutronics_reaction_rates.size(); ++nuclide) { - unsigned int num_fission = std::floor(recoil_rate_scaling * dt * vol * - (*_partial_neutronics_reaction_rates[nuclide]) / (*_averaged_number_densities[nuclide] * averaged_data._site_volume) - * averaged_data._elements[nuclide] + getRandomReal()); + unsigned int num_fission = + std::floor(recoil_rate_scaling * dt * vol * (*_partial_neutronics_reaction_rates[nuclide]) / + (*_averaged_number_densities[nuclide] * averaged_data._site_volume) * + averaged_data._elements[nuclide] + + getRandomReal()); for (unsigned i = 0; i < num_fission; ++i) { @@ -59,8 +72,8 @@ PKAFissionFragmentNeutronics::appendPKAs(std::vector & 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(pka_parameters, ion[0]._Z, ion[0]._m); + ion[1]._tag = ionTag(pka_parameters, ion[1]._Z, ion[1]._m); // set location of the fission event setPosition(ion[0]); diff --git a/src/userobjects/PKAFixedPointGenerator.C b/src/userobjects/PKAFixedPointGenerator.C index 0046282c..c7d140ab 100644 --- a/src/userobjects/PKAFixedPointGenerator.C +++ b/src/userobjects/PKAFixedPointGenerator.C @@ -11,12 +11,15 @@ registerMooseObject("MagpieApp", PKAFixedPointGenerator); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addClassDescription("This PKAGenerator starts particle from a fixed point in a random direction (isotropically)."); - params.addParam("num_pkas", 1000, "The number of PKAs to be started from this position"); + params.addClassDescription("This PKAGenerator starts particle from a fixed point in a random " + "direction (isotropically)."); + params.addParam( + "num_pkas", 1000, "The number of PKAs to be started from this position"); params.addRequiredParam("point", "The point from which the PKAs are started"); params.addRequiredParam("Z", "PKA nuclear charge"); params.addRequiredParam("m", "PKA mass in amu"); @@ -24,8 +27,8 @@ InputParameters validParams() return params; } -PKAFixedPointGenerator::PKAFixedPointGenerator(const InputParameters & parameters) : - PKAGeneratorBase(parameters), +PKAFixedPointGenerator::PKAFixedPointGenerator(const InputParameters & parameters) + : PKAGeneratorBase(parameters), _num_pka(getParam("num_pkas")), _point(getParam("point")), _Z(getParam("Z")), @@ -38,14 +41,16 @@ PKAFixedPointGenerator::PKAFixedPointGenerator(const InputParameters & parameter } void -PKAFixedPointGenerator::appendPKAs(std::vector & ion_list, Real /*dt*/, Real /*vol*/, Real recoil_rate_scaling, const MyTRIMRasterizer::AveragedData & averaged_data) const +PKAFixedPointGenerator::appendPKAs(std::vector & ion_list, + const MyTRIMRasterizer::PKAParameters & pka_parameters, + const MyTRIMRasterizer::AveragedData &) const { if (_current_elem->id() != _elem_id) return; unsigned int num_pka = _num_pka; - if (recoil_rate_scaling != 1) - num_pka = std::floor(recoil_rate_scaling * _num_pka + getRandomReal()); + if (pka_parameters._recoil_rate_scaling != 1) + num_pka = std::floor(pka_parameters._recoil_rate_scaling * _num_pka + getRandomReal()); for (unsigned i = 0; i < _num_pka; ++i) { @@ -59,7 +64,8 @@ PKAFixedPointGenerator::appendPKAs(std::vector & 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(pka_parameters, pka._Z, pka._m); + ; // set stopping criteria pka.setEf(); diff --git a/src/userobjects/PKAGeneratorAlphaDecay.C b/src/userobjects/PKAGeneratorAlphaDecay.C index a182a5c3..c0fcc9db 100644 --- a/src/userobjects/PKAGeneratorAlphaDecay.C +++ b/src/userobjects/PKAGeneratorAlphaDecay.C @@ -12,20 +12,24 @@ registerMooseObject("MagpieApp", PKAGeneratorAlphaDecay); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addClassDescription("A PKAGenerator for starting alpha particles from decay\nDecay data is retrieved from file data/alpha_decay/alpha_decay.txt"); - params.addRequiredParam>("ZAID", "The Z/A ids corresponding to the var arguments in the rasterizer [1e4 * A + 10 * Z + state]"); + params.addClassDescription("A PKAGenerator for starting alpha particles from decay\nDecay data " + "is retrieved from file data/alpha_decay/alpha_decay.txt"); + params.addRequiredParam>("ZAID", + "The Z/A ids corresponding to the var " + "arguments in the rasterizer [1e4 * A + 10 * " + "Z + state]"); MooseEnum timeUnit("second=0 millisecond=1 microsecond=2", "second"); params.addParam("time_unit", timeUnit, "Unit of time used in this model"); return params; } -PKAGeneratorAlphaDecay::PKAGeneratorAlphaDecay(const InputParameters & parameters) : - PKAGeneratorBase(parameters), - _zaids(getParam>("ZAID")) +PKAGeneratorAlphaDecay::PKAGeneratorAlphaDecay(const InputParameters & parameters) + : PKAGeneratorBase(parameters), _zaids(getParam>("ZAID")) { switch (getParam("time_unit")) { @@ -97,25 +101,32 @@ PKAGeneratorAlphaDecay::readAlphaData() } } - // in debug mode we want some info on what is stored in _decay_data_sets - #if DEBUG - _console << "\nAlpha decay read from file\n"; - for (auto & d : _decay_data_sets) - { - _console << "ZAID: " << d.first << " has " << d.second.size() << " decay channels:\n"; - for (unsigned int j = 0; j < d.second.size(); ++j) - _console << "lambda = " << d.second[j]._decay_constants << " " - << "energy = " << d.second[j]._alpha_energies << " " - << "intensity = " << d.second[j]._intensities << "\n"; - } - _console << std::endl; - #endif +// in debug mode we want some info on what is stored in _decay_data_sets +#if DEBUG + _console << "\nAlpha decay read from file\n"; + for (auto & d : _decay_data_sets) + { + _console << "ZAID: " << d.first << " has " << d.second.size() << " decay channels:\n"; + for (unsigned int j = 0; j < d.second.size(); ++j) + _console << "lambda = " << d.second[j]._decay_constants << " " + << "energy = " << d.second[j]._alpha_energies << " " + << "intensity = " << d.second[j]._intensities << "\n"; + } + _console << std::endl; +#endif } void -PKAGeneratorAlphaDecay::appendPKAs(std::vector & ion_list, Real dt, Real vol, Real recoil_rate_scaling, const MyTRIMRasterizer::AveragedData & averaged_data) const +PKAGeneratorAlphaDecay::appendPKAs(std::vector & ion_list, + const MyTRIMRasterizer::PKAParameters & pka_parameters, + const MyTRIMRasterizer::AveragedData & averaged_data) const { - mooseAssert(dt >= 0, "Passed a negative time window into PKAFissionFragmentNeutronics::appendPKAs"); + const auto dt = pka_parameters._dt; + const auto vol = pka_parameters._volume; + const auto recoil_rate_scaling = pka_parameters._recoil_rate_scaling; + + mooseAssert(dt >= 0, + "Passed a negative time window into PKAFissionFragmentNeutronics::appendPKAs"); mooseAssert(vol >= 0, "Passed a negative volume into PKAFissionFragmentNeutronics::appendPKAs"); if (averaged_data._elements.size() != _zaids.size()) @@ -133,8 +144,10 @@ PKAGeneratorAlphaDecay::appendPKAs(std::vector & ion_list, R auto & decay = it->second; for (unsigned int l = 0; l < decay.size(); ++l) { - unsigned int num_decay = std::floor(recoil_rate_scaling * vol * decay[l]._intensities * averaged_data._elements[nuclide] - * (1.0 - std::exp(-decay[l]._decay_constants * dt)) + getRandomReal()); + unsigned int num_decay = std::floor(recoil_rate_scaling * vol * decay[l]._intensities * + averaged_data._elements[nuclide] * + (1.0 - std::exp(-decay[l]._decay_constants * dt)) + + getRandomReal()); for (unsigned i = 0; i < num_decay; ++i) { @@ -155,8 +168,8 @@ PKAGeneratorAlphaDecay::appendPKAs(std::vector & 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(pka_parameters, ion[0]._Z, ion[0]._m); + ion[1]._tag = ionTag(pka_parameters, ion[1]._Z, ion[1]._m); // set the energies; use momentum conservation ion[0]._E = decay[l]._alpha_energies; diff --git a/src/userobjects/PKAGeneratorBase.C b/src/userobjects/PKAGeneratorBase.C index 0841e27d..8887cb3d 100644 --- a/src/userobjects/PKAGeneratorBase.C +++ b/src/userobjects/PKAGeneratorBase.C @@ -10,16 +10,18 @@ #include "MagpieUtils.h" #include -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - 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,23 +33,31 @@ PKAGeneratorBase::setPosition(MyTRIM_NS::IonBase & ion) const } int -PKAGeneratorBase::ionTag(const std::vector & rasterizer_Z, const std::vector & rasterizer_m, Real Z, Real m) const +PKAGeneratorBase::ionTag(const MyTRIMRasterizer::PKAParameters & pka_parameters, + 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; - return -1; + // this function relies on the exact representation of whole numbers in IEEE floating point + // numbers up to a reasonable upper limit [Z < m < 300] + + // element not found in rasterizer table + if (pka_parameters._num_Z[Z] == 0) + return -1; + + if (pka_parameters._num_Z[Z] == 1) + return pka_parameters._single_Z_index[Z]; + + const auto & mZ = pka_parameters._mass_charge_pair; + // start the search at the firsat matching Z + for (auto i = pka_parameters._single_Z_index[Z]; i < mZ.size(); ++i) + if (mZ[i].second == Z && mZ[i].first == m) + return i; + + mooseError("Unexpectedly did not find the element with Z=", + Z, + " m=", + m, + ". Inconsistency in _numZ and _single_Z_index."); } void diff --git a/src/userobjects/PKAGeneratorRecoil.C b/src/userobjects/PKAGeneratorRecoil.C index 9c1a2144..2310b903 100644 --- a/src/userobjects/PKAGeneratorRecoil.C +++ b/src/userobjects/PKAGeneratorRecoil.C @@ -12,31 +12,42 @@ registerMooseObject("MagpieApp", PKAGeneratorRecoil); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addClassDescription("PKA recoil generator user object.\n Takes pdf and samples PKAs due to recoil reaction."); + params.addClassDescription( + "PKA recoil generator user object.\n Takes pdf and samples PKAs due to recoil reaction."); return params; } -PKAGeneratorRecoil::PKAGeneratorRecoil(const InputParameters & parameters): - PKAGeneratorNeutronicsBase(parameters) +PKAGeneratorRecoil::PKAGeneratorRecoil(const InputParameters & parameters) + : PKAGeneratorNeutronicsBase(parameters) { } void -PKAGeneratorRecoil::setPDF(const std::vector & ZAID, const std::vector & energies, const MultiIndex & probabilities) +PKAGeneratorRecoil::setPDF(const std::vector & ZAID, + const std::vector & energies, + const MultiIndex & probabilities) { _pdf = DiscretePKAPDF(ZAID, energies, probabilities); - #if DEBUG - _console << "\nPKAGeneratorRecoil object received the following pdf from DiscretePKAPDF object:\n" << _pdf; - #endif +#if DEBUG + _console << "\nPKAGeneratorRecoil object received the following pdf from DiscretePKAPDF object:\n" + << _pdf; +#endif } void -PKAGeneratorRecoil::appendPKAs(std::vector & ion_list, Real dt, Real vol, Real recoil_rate_scaling, const MyTRIMRasterizer::AveragedData & averaged_data) const +PKAGeneratorRecoil::appendPKAs(std::vector & ion_list, + const MyTRIMRasterizer::PKAParameters & pka_parameters, + const MyTRIMRasterizer::AveragedData & averaged_data) const { + const auto dt = pka_parameters._dt; + const auto vol = pka_parameters._volume; + const auto recoil_rate_scaling = pka_parameters._recoil_rate_scaling; + mooseAssert(dt >= 0, "Passed a negative time window into PKAGeneratorRecoil::appendPKAs"); mooseAssert(vol >= 0, "Passed a negative volume into PKAGeneratorRecoil::appendPKAs"); @@ -45,9 +56,10 @@ PKAGeneratorRecoil::appendPKAs(std::vector & ion_list, Real for (unsigned int nuclide = 0; nuclide < _partial_neutronics_reaction_rates.size(); ++nuclide) { - unsigned int num_recoils = std::floor(recoil_rate_scaling * dt * vol * - (*_partial_neutronics_reaction_rates[nuclide]) / (*_averaged_number_densities[nuclide]) - * averaged_data._elements[nuclide] + getRandomReal()); + unsigned int num_recoils = + std::floor(recoil_rate_scaling * dt * vol * (*_partial_neutronics_reaction_rates[nuclide]) / + (*_averaged_number_densities[nuclide]) * averaged_data._elements[nuclide] + + getRandomReal()); for (unsigned i = 0; i < num_recoils; ++i) { diff --git a/src/userobjects/PKAGun.C b/src/userobjects/PKAGun.C index ed2ea811..2f0ad0c1 100644 --- a/src/userobjects/PKAGun.C +++ b/src/userobjects/PKAGun.C @@ -10,18 +10,19 @@ registerMooseObject("MagpieApp", PKAGun); -template<> -InputParameters validParams() +template <> +InputParameters +validParams() { InputParameters params = validParams(); - params.addClassDescription("This PKAGenerator starts particle from a fixed point in a fixed direction."); + params.addClassDescription( + "This PKAGenerator starts particle from a fixed point in a fixed direction."); params.addRequiredParam("direction", "The fixed direction the PKAs move along"); return params; } -PKAGun::PKAGun(const InputParameters & parameters) : - PKAFixedPointGenerator(parameters), - _direction(getParam("direction")) +PKAGun::PKAGun(const InputParameters & parameters) + : PKAFixedPointGenerator(parameters), _direction(getParam("direction")) { } From 3935d4f24bc3eaa441956ff46287a15cb72f16e9 Mon Sep 17 00:00:00 2001 From: Daniel Schwen Date: Mon, 13 Aug 2018 10:36:09 -0600 Subject: [PATCH 2/3] Redo Z indexing (#333) --- include/userobjects/MyTRIMRasterizer.h | 8 +++----- src/userobjects/MyTRIMRasterizer.C | 22 ++++++++++++++++------ src/userobjects/PKAGeneratorBase.C | 8 ++++---- 3 files changed, 23 insertions(+), 15 deletions(-) diff --git a/include/userobjects/MyTRIMRasterizer.h b/include/userobjects/MyTRIMRasterizer.h index 4695be6b..88a51ec3 100644 --- a/include/userobjects/MyTRIMRasterizer.h +++ b/include/userobjects/MyTRIMRasterizer.h @@ -69,11 +69,9 @@ class MyTRIMRasterizer : public ElementUserObject /// masses charges (m_i, Z_i) of the matrix elements std::vector> _mass_charge_pair; - /// how many isotopes to we have for each element? (only support up to Z=119!) - std::array _num_Z; - - /// if only one element with a specific Z is present point to its rasterizer index here for fast lookup - std::array _single_Z_index; + /// how many isotopes to we have for each element and which index contains the + /// first Z match? (only support up to Z=119!) + std::array, 120> _index_Z; /// time interval over which the PKAs are added Real _dt; diff --git a/src/userobjects/MyTRIMRasterizer.C b/src/userobjects/MyTRIMRasterizer.C index f831735a..0877a813 100644 --- a/src/userobjects/MyTRIMRasterizer.C +++ b/src/userobjects/MyTRIMRasterizer.C @@ -15,6 +15,8 @@ #include "libmesh/quadrature.h" #include "libmesh/parallel_algebra.h" +#include + // custom data load and data store methods for a struct with an std::vector member template <> inline void @@ -163,9 +165,17 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) if (trim_Z.size() != _nvars) mooseError("Parameter 'Z' must have as many components as coupled variables."); + // error check masses and charges 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); + if (trim_Z[i] >= _pka_parameters._index_Z.size()) + mooseError("Value of Z is too large. Maximum Z supported is ", + _pka_parameters._index_Z.size() - 1, + " but one element has Z=", + i); + } for (unsigned int i = 0; i < _nvars; ++i) { @@ -240,8 +250,8 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) } // setup invariant PKA generation parameters - for (auto & nZ : _pka_parameters._num_Z) - nZ = 0; + for (auto & nZ : _pka_parameters._index_Z) + nZ = std::make_pair(0, 0); _pka_parameters._mass_charge_pair.resize(_nvars); _pka_parameters._recoil_rate_scaling = _trim_parameters.recoil_rate_scaling; for (unsigned int i = 0; i < _nvars; ++i) @@ -253,12 +263,12 @@ MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters) std::make_pair(_trim_parameters.element_prototypes[i]._m, Z); // increase the count of elements with the same Z - auto & num_Z = _pka_parameters._num_Z[Z]; - num_Z++; + auto & index_Z = _pka_parameters._index_Z[Z]; + index_Z.first++; // only set this to the first index (important for ionTag()) - if (num_Z == 1) - _pka_parameters._single_Z_index[Z] = i; + if (index_Z.first == 1) + index_Z.second = i; } } diff --git a/src/userobjects/PKAGeneratorBase.C b/src/userobjects/PKAGeneratorBase.C index 8887cb3d..f43e3ed2 100644 --- a/src/userobjects/PKAGeneratorBase.C +++ b/src/userobjects/PKAGeneratorBase.C @@ -41,15 +41,15 @@ PKAGeneratorBase::ionTag(const MyTRIMRasterizer::PKAParameters & pka_parameters, // numbers up to a reasonable upper limit [Z < m < 300] // element not found in rasterizer table - if (pka_parameters._num_Z[Z] == 0) + if (pka_parameters._index_Z[Z].first == 0) return -1; - if (pka_parameters._num_Z[Z] == 1) - return pka_parameters._single_Z_index[Z]; + if (pka_parameters._index_Z[Z].first == 1) + return pka_parameters._index_Z[Z].second; const auto & mZ = pka_parameters._mass_charge_pair; // start the search at the firsat matching Z - for (auto i = pka_parameters._single_Z_index[Z]; i < mZ.size(); ++i) + for (auto i = pka_parameters._index_Z[Z].second; i < mZ.size(); ++i) if (mZ[i].second == Z && mZ[i].first == m) return i; From d791a47b8e33d11b98d43afb7999a178cd8cd8e2 Mon Sep 17 00:00:00 2001 From: Daniel Schwen Date: Mon, 13 Aug 2018 12:14:37 -0600 Subject: [PATCH 3/3] Fix failed match bug (#333) --- src/userobjects/PKAGeneratorBase.C | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/userobjects/PKAGeneratorBase.C b/src/userobjects/PKAGeneratorBase.C index f43e3ed2..bbf99f74 100644 --- a/src/userobjects/PKAGeneratorBase.C +++ b/src/userobjects/PKAGeneratorBase.C @@ -53,11 +53,8 @@ PKAGeneratorBase::ionTag(const MyTRIMRasterizer::PKAParameters & pka_parameters, if (mZ[i].second == Z && mZ[i].first == m) return i; - mooseError("Unexpectedly did not find the element with Z=", - Z, - " m=", - m, - ". Inconsistency in _numZ and _single_Z_index."); + // no matching mass (isotope) found for the given Z + return -1; } void