diff --git a/doc/content/bib/magpie.bib b/doc/content/bib/magpie.bib new file mode 100644 index 00000000..ef598ba8 --- /dev/null +++ b/doc/content/bib/magpie.bib @@ -0,0 +1,8 @@ + +@inproceedings{SchunertTREATHeatSource, + author={S. Schunert and others}, + title={{Heat Source Characterization In A TREAT Fuel Particle Using Coupled Neutronics Binary Collision Monte-Carlo Calculations}}, + booktitle={{M\&C 2017 International Conference on Mathematics \& Computational Methods Applied to Nuclear Science \& Engineering}}, + month={April}, + year={2017} +} diff --git a/doc/content/documentation/systems/UserObjects/PKAFissionFragmentNeutronics.md b/doc/content/documentation/systems/UserObjects/PKAFissionFragmentNeutronics.md index 388394a6..477c8221 100644 --- a/doc/content/documentation/systems/UserObjects/PKAFissionFragmentNeutronics.md +++ b/doc/content/documentation/systems/UserObjects/PKAFissionFragmentNeutronics.md @@ -1,12 +1,54 @@ - - # PKAFissionFragmentNeutronics -!alert construction title=Undocumented Class -The PKAFissionFragmentNeutronics has not been documented, if you would like to contribute to MOOSE by -writing documentation, please see [/generate.md]. The content contained on this page explains -the typical documentation associated with a MooseObject; however, what is contained is ultimately -determined by what is necessary to make the documentation clear for users. +`PKAFissionFragmentNeutronics` generates primary knock-on atoms (PKA) originating +from fission reactions. The difference to PKAFissionFragmentEmpirical is that +it uses isotopic fission rates computed by neutronics calculations and ENDF data +for sampling fission product species. + +The partial fission rate for isotope $i$ is denoted by $F_i(\vec{r}, \vec{\rho})$. +It is computed by: + +\begin{equation} +F_i(\vec{r}, \vec{\rho}) = \sum\limits_{g=1}^G N_i(\vec{r},\vec{\rho})\sigma_{f,g,i} \phi_g(\vec{r}), +\end{equation} + +where $\vec{r}$ is the location in the macroscopic, neutronics domain and $\vec{\rho}$ +is the location in the microscopic domain. These two locations are separated because their +scale is significantly separated by 3-4 orders of magnitude. Changes with $\vec{r}$ are +understood to be smooth changes of the average composition of the microscopic domain (orders of centimeters), while changes with $\vec{\rho}$ captures compositional changes between grains on the microscopic domain, i.e. on the orders of micro-meters. The notation of separating the spatial dependence into two scales follows homogenization theory. $N_i$ is the number density of isotope $i$, $\sigma_{f,g,i}$ is the microscopic fission cross section of isotope $i$ in group $g$, and $\phi_g$ is the scalar flux in group $g$. + +Denoting the microscopic domain as $\Omega_m$ located at $\vec{r}_m$, we can define a slowly varying average of the number densities: + +\begin{equation} + N_i(\vec{r}) = \frac{1}{\Omega_m} \int_{\Omega_m} N_i(\vec{r},\vec{\rho}) d\vec{\rho}. +\end{equation} + +$N_i(\vec{r})$ is the number density provided to neutronics calculations. The nuclide fission rates computed by the neutronics calculation is consequently the slowly varying average +given by: + +\begin{equation} +F_i (\vec{r}) = \sum\limits_{g=1}^G N_i(\vec{r})\sigma_{f,g,i}. +\end{equation} + +The `PKAFissionFragmentNeutronics` accepts the values of $F_i (\vec{r}_m)$ as the `partial_reaction_rates` parameter. In addition, it accepts the values of $N_i(\vec{r}_m)$ +as the `averaged_number_densities` parameter. + +The fission rate density at a location in the microscopic domain is computed by: + +\begin{equation} + F_i (\vec{r}_m, \vec{\rho}) = \frac{N_i(\vec{r}_m, \vec{\rho})}{N_i(\vec{r}_m)} F_i(\vec{r}_m), +\end{equation} + +where $N_i(\vec{r}_m, \vec{\rho}) $ are the number densities provided by the rasterizer. + +The expected number of fissions in a mesh element with index $j$ and volume $V_j$ is given by: + +\begin{equation} +C_i = \Delta t \int_{V_j} F_i (\vec{r}_m, \vec{\rho}) d\vec{\rho}. +\end{equation} + +Non-integer results are rounded up with a probability of the $C_i - \text{int}(C_i)$; otherwise +rounded down. Two PKAs are created from each fission event. The algorithm for sampling their type, energy, and direction of motion is described in [cite:SchunertTREATHeatSource]. !syntax description /UserObjects/PKAFissionFragmentNeutronics diff --git a/doc/content/index.md b/doc/content/index.md index 74bebb96..6cdbb41d 100644 --- a/doc/content/index.md +++ b/doc/content/index.md @@ -7,7 +7,7 @@

Mesoscale-Atomistics Glue Program for Integrated Execution

Welcome to the Magpie project. Magpie is an Idaho National Laboratory tool for coupling -atomistic and Monte Carlo methods to the finite element MOOSE framework. Iy currently couples -the MyTRIM binary collision Montecarlo code and teh Sandia SPPARKS Monte Carlo toolkit. +atomistic and Monte Carlo methods to the finite element MOOSE framework. It currently couples +the MyTRIM binary collision Monte-Carlo code and the Sandia SPPARKS Monte Carlo toolkit. diff --git a/doc/hidden.yml b/doc/hidden.yml index 27082f04..3e91eae6 100644 --- a/doc/hidden.yml +++ b/doc/hidden.yml @@ -21,7 +21,6 @@ - /UserObjects/NeutronicsSpectrumSamplerFission - /UserObjects/PKAConstant - /UserObjects/PKAFissionFragmentEmpirical -- /UserObjects/PKAFissionFragmentNeutronics - /UserObjects/PKAFixedPointGenerator - /UserObjects/PKAFunction - /UserObjects/PKAGeneratorAlphaDecay diff --git a/include/userobjects/PKAFissionFragmentNeutronics.h b/include/userobjects/PKAFissionFragmentNeutronics.h index 6536398a..8d77d9c8 100644 --- a/include/userobjects/PKAFissionFragmentNeutronics.h +++ b/include/userobjects/PKAFissionFragmentNeutronics.h @@ -34,17 +34,6 @@ class PKAFissionFragmentNeutronics : public PKAGeneratorNeutronicsBase * will be sampled to determine the initial state of the ions passed to MyTRIM */ DiscreteFissionPKAPDF _pdf; - - /** - * _partial_fission_rates: the fission rate per nuclide is: - * f_i = N_i * (sum_g sigma_{f,g,i} phi_g), - * where i is the nuclide ID and g is the energy group. - * N_i is provided as variable in the Rasterizer, but the (sum_g sigma_{f,g,i} phi_g) must - * be provided as partial fission rate here - * Note: the fission rate densities are in units of [ fissions / ( mesh-length-unit^3 sec ) ] - */ - std::vector _partial_fission_rates; - std::vector _stored_pps; }; #endif //PKAFISSIONFRAGMENTNEUTRONICS_H diff --git a/include/userobjects/PKAGeneratorNeutronicsBase.h b/include/userobjects/PKAGeneratorNeutronicsBase.h index 418bdbc2..48bb678b 100644 --- a/include/userobjects/PKAGeneratorNeutronicsBase.h +++ b/include/userobjects/PKAGeneratorNeutronicsBase.h @@ -31,6 +31,20 @@ class PKAGeneratorNeutronicsBase : public PKAGeneratorBase /// helper function to set pdf based on neutronics information from Mammoth virtual void setPDF(const std::vector & ZAID, const std::vector & energies, const MultiIndex & probabilities) = 0; + +protected: + /** + * The partial reaction rates are reaction rates of a certain type from the neutronics calculation + * for a single nculide, e.g. fission rate of U-235. NOTE: it multiplies the smoothly varying _average_ + * of the U-235 number density already + * + * In Magpie, the rate of PKA creation is obtained by multiplying with the _local_ number density. To obtain a + * proper PKA creation rate, we need to divide by the smooth average. This is what _averaged_number_densities is for. + */ + std::vector _partial_neutronics_reaction_rates; + std::vector _stored_reaction_rates; + std::vector _averaged_number_densities; + std::vector _stored_densities; }; #endif // PKAGENERATORNEUTRONICSBASE_H diff --git a/include/userobjects/PKAGeneratorRecoil.h b/include/userobjects/PKAGeneratorRecoil.h index 445ab21b..3f3e9e93 100644 --- a/include/userobjects/PKAGeneratorRecoil.h +++ b/include/userobjects/PKAGeneratorRecoil.h @@ -33,13 +33,6 @@ class PKAGeneratorRecoil : public PKAGeneratorNeutronicsBase * are sampled */ DiscretePKAPDF _pdf; - - /** - * The partial reaction rates [# reactions / time / volume] of the particular - * recoil reaction for each nuclide - */ - std::vector _partial_recoil_rates; - std::vector _stored_pps; }; #endif //PKAGeneratorRecoil_H diff --git a/src/userobjects/PKAFissionFragmentNeutronics.C b/src/userobjects/PKAFissionFragmentNeutronics.C index e6814d95..d5805869 100644 --- a/src/userobjects/PKAFissionFragmentNeutronics.C +++ b/src/userobjects/PKAFissionFragmentNeutronics.C @@ -17,8 +17,6 @@ template<> InputParameters validParams() { InputParameters params = validParams(); - params.addParam>("partial_fission_rates", "Partial fission rates per unit volume [sum_g sigma_{f,g,i} * phi_g]. " - "Provide number density as variable in rasterizer!"); params.addClassDescription("PKA generator (fission) user object.\n Takes pdf and samples PKAs due to fission."); return params; } @@ -26,32 +24,6 @@ InputParameters validParams() PKAFissionFragmentNeutronics::PKAFissionFragmentNeutronics(const InputParameters & parameters): PKAGeneratorNeutronicsBase(parameters) { - if (isParamValid("partial_fission_rates")) - { - std::vector names = getParam>("partial_fission_rates"); - _partial_fission_rates.resize(names.size()); - _stored_pps.resize(names.size()); - for (unsigned int j = 0; j < names.size(); ++j) - if (_fe_problem.hasPostprocessor(names[j])) - _partial_fission_rates[j] = &getPostprocessorValueByName(names[j]); - else - { - Real real_value = -std::numeric_limits::max(); - std::istringstream ss(names[j]); - - if (ss >> real_value && ss.eof()) - _stored_pps[j] = real_value; - else - mooseError("Illegal entry in partial_fission_rates: ", names[j]); - - _partial_fission_rates[j] = &_stored_pps[j]; - } - } - else - { - _stored_pps = {1.0e-8}; - _partial_fission_rates = {&_stored_pps[0]}; - } } void @@ -66,12 +38,14 @@ PKAFissionFragmentNeutronics::appendPKAs(std::vector & ion_l 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_fission_rates.size()) - mooseError("Size of averaged_data and partial_fission_rates must be equal"); + if (averaged_data._elements.size() != _partial_neutronics_reaction_rates.size()) + mooseError("Size of averaged_data and partial_reaction_rates must be equal"); - for (unsigned int nuclide = 0; nuclide < _partial_fission_rates.size(); ++nuclide) + for (unsigned int nuclide = 0; nuclide < _partial_neutronics_reaction_rates.size(); ++nuclide) { - unsigned int num_fission = std::floor(recoil_rate_scaling * dt * vol * (*_partial_fission_rates[nuclide]) * 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._elements[nuclide] + getRandomReal()); for (unsigned i = 0; i < num_fission; ++i) { diff --git a/src/userobjects/PKAGeneratorNeutronicsBase.C b/src/userobjects/PKAGeneratorNeutronicsBase.C index 80521f80..8fa54de7 100644 --- a/src/userobjects/PKAGeneratorNeutronicsBase.C +++ b/src/userobjects/PKAGeneratorNeutronicsBase.C @@ -15,6 +15,9 @@ template<> InputParameters validParams() { InputParameters params = validParams(); + params.addParam>("partial_reaction_rates", "Partial neutronic reaction rates per unit volume [sum_g xs_{r,g,i} * phi_g], " + "r: reaction type, g: energy group, i: nuclide id.Provide number density as variable in rasterizer!"); + params.addParam>("averaged_number_densities", "The number density of the species averaged over the domain."); params.addClassDescription("PKA generator (neutronics) user object base class.\n Takes pdf and samples PKAs due to various interactions."); return params; } @@ -22,4 +25,65 @@ InputParameters validParams() PKAGeneratorNeutronicsBase::PKAGeneratorNeutronicsBase(const InputParameters & parameters) : PKAGeneratorBase(parameters) { + if (isParamValid("partial_reaction_rates")) + { + std::vector names = getParam>("partial_reaction_rates"); + _partial_neutronics_reaction_rates.resize(names.size()); + _stored_reaction_rates.resize(names.size()); + for (unsigned int j = 0; j < names.size(); ++j) + if (_fe_problem.hasPostprocessor(names[j])) + _partial_neutronics_reaction_rates[j] = &getPostprocessorValueByName(names[j]); + else + { + Real real_value = -std::numeric_limits::max(); + std::istringstream ss(names[j]); + + if (ss >> real_value && ss.eof()) + _stored_reaction_rates[j] = real_value; + else + mooseError("Illegal entry in _partial_neutronics_reaction_rates: ", names[j]); + + _partial_neutronics_reaction_rates[j] = &_stored_reaction_rates[j]; + } + } + else + { + _stored_reaction_rates = {1.0e-8}; + _partial_neutronics_reaction_rates = {&_stored_reaction_rates[0]}; + } + + if (isParamValid("averaged_number_densities")) + { + std::vector names = getParam>("averaged_number_densities"); + _averaged_number_densities.resize(names.size()); + _stored_densities.resize(names.size()); + for (unsigned int j = 0; j < names.size(); ++j) + if (_fe_problem.hasPostprocessor(names[j])) + _averaged_number_densities[j] = &getPostprocessorValueByName(names[j]); + else + { + Real real_value = -std::numeric_limits::max(); + std::istringstream ss(names[j]); + + if (ss >> real_value && ss.eof()) + _stored_densities[j] = real_value; + else + mooseError("Illegal entry in _partial_neutronics_reaction_rates: ", names[j]); + + _averaged_number_densities[j] = &_stored_reaction_rates[j]; + } + } + else + { + _averaged_number_densities.resize(_stored_reaction_rates.size()); + _stored_densities.resize(_stored_reaction_rates.size()); + for (unsigned int j = 0; j < _stored_densities.size(); ++j) + { + _stored_densities[j] = 1; + _averaged_number_densities[j] = &_stored_densities[j]; + } + } + + if (_averaged_number_densities.size() != _partial_neutronics_reaction_rates.size()) + mooseError("partial_reaction_rates and averaged_number_densities must have the same number of entries."); } diff --git a/src/userobjects/PKAGeneratorRecoil.C b/src/userobjects/PKAGeneratorRecoil.C index fd791c40..9c1a2144 100644 --- a/src/userobjects/PKAGeneratorRecoil.C +++ b/src/userobjects/PKAGeneratorRecoil.C @@ -16,8 +16,6 @@ template<> InputParameters validParams() { InputParameters params = validParams(); - params.addParam>("partial_recoil_rates", "Partial recoil rates per unit volume [recoil reaction rate per nuclide]. " - "Provide number density as variable in rasterizer!"); params.addClassDescription("PKA recoil generator user object.\n Takes pdf and samples PKAs due to recoil reaction."); return params; } @@ -25,32 +23,6 @@ InputParameters validParams() PKAGeneratorRecoil::PKAGeneratorRecoil(const InputParameters & parameters): PKAGeneratorNeutronicsBase(parameters) { - if (isParamValid("partial_recoil_rates")) - { - std::vector names = getParam>("partial_recoil_rates"); - _partial_recoil_rates.resize(names.size()); - _stored_pps.resize(names.size()); - for (unsigned int j = 0; j < names.size(); ++j) - if (_fe_problem.hasPostprocessor(names[j])) - _partial_recoil_rates[j] = &getPostprocessorValueByName(names[j]); - else - { - Real real_value = -std::numeric_limits::max(); - std::istringstream ss(names[j]); - - if (ss >> real_value && ss.eof()) - _stored_pps[j] = real_value; - else - mooseError("Illegal entry in partial_recoil_rates: ", names[j]); - - _partial_recoil_rates[j] = &_stored_pps[j]; - } - } - else - { - _stored_pps = {1.0e-8}; - _partial_recoil_rates = {&_stored_pps[0]}; - } } void @@ -68,12 +40,14 @@ PKAGeneratorRecoil::appendPKAs(std::vector & ion_list, Real mooseAssert(dt >= 0, "Passed a negative time window into PKAGeneratorRecoil::appendPKAs"); mooseAssert(vol >= 0, "Passed a negative volume into PKAGeneratorRecoil::appendPKAs"); - if (averaged_data._elements.size() != _partial_recoil_rates.size()) - mooseError("Size of averaged_data and partial_recoil_rates must be equal"); + if (averaged_data._elements.size() != _partial_neutronics_reaction_rates.size()) + mooseError("Size of averaged_data and partial_reaction_rates must be equal"); - for (unsigned int nuclide = 0; nuclide < _partial_recoil_rates.size(); ++nuclide) + for (unsigned int nuclide = 0; nuclide < _partial_neutronics_reaction_rates.size(); ++nuclide) { - unsigned int num_recoils = std::floor(recoil_rate_scaling * dt * vol * (*_partial_recoil_rates[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/tests/radiation_damage/coupled_fission_mockup/damage_sub.i b/tests/radiation_damage/coupled_fission_mockup/damage_sub.i index caa11f37..98da393b 100644 --- a/tests/radiation_damage/coupled_fission_mockup/damage_sub.i +++ b/tests/radiation_damage/coupled_fission_mockup/damage_sub.i @@ -137,7 +137,7 @@ [./neutronics_fission_generator] type = PKAFissionFragmentNeutronics relative_density = 1 - partial_fission_rates = '2.0e-11 0 0' + partial_reaction_rates = '2.0e-11 0 0' [../] [./rasterizer] type = MyTRIMRasterizer diff --git a/tests/radiation_damage/coupled_fission_mockup/damage_sub_3D.i b/tests/radiation_damage/coupled_fission_mockup/damage_sub_3D.i index 86492849..40cc18bc 100644 --- a/tests/radiation_damage/coupled_fission_mockup/damage_sub_3D.i +++ b/tests/radiation_damage/coupled_fission_mockup/damage_sub_3D.i @@ -138,7 +138,7 @@ [./neutronics_fission_generator] type = PKAFissionFragmentNeutronics relative_density = 1 - partial_fission_rates = '1.0e-14 0 0' + partial_reaction_rates = '1.0e-14 0 0' [../] [./rasterizer] type = MyTRIMRasterizer