From ce02d5fefd5d40755d8105eb771427d75dc8ab9e Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Mon, 13 Jan 2020 14:24:57 -0700 Subject: [PATCH] Implement dpa userobject; inherit from NeutronDamageInterface to allow different damage functions to be plugged in (#410) --- include/postprocessors/DPAPostprocessor.h | 34 ++ include/userobjects/DPAUserObjectBase.h | 101 +++++ include/userobjects/FunctionDPAUserObject.h | 44 ++ .../userobjects/ParkinCoulterDPAUserObject.h | 38 ++ ...CoulterBase.h => ParkinCoulterInterface.h} | 16 +- include/userobjects/PolyatomicRecoil.h | 8 +- src/postprocessors/DPAPostprocessor.C | 37 ++ src/userobjects/DPAUserObjectBase.C | 406 ++++++++++++++++++ src/userobjects/FunctionDPAUserObject.C | 151 +++++++ src/userobjects/ParkinCoulterDPAUserObject.C | 94 ++++ ...CoulterBase.C => ParkinCoulterInterface.C} | 35 +- src/userobjects/PolyatomicRecoil.C | 5 +- 12 files changed, 938 insertions(+), 31 deletions(-) create mode 100644 include/postprocessors/DPAPostprocessor.h create mode 100644 include/userobjects/DPAUserObjectBase.h create mode 100644 include/userobjects/FunctionDPAUserObject.h create mode 100644 include/userobjects/ParkinCoulterDPAUserObject.h rename include/userobjects/{ParkinCoulterBase.h => ParkinCoulterInterface.h} (81%) create mode 100644 src/postprocessors/DPAPostprocessor.C create mode 100644 src/userobjects/DPAUserObjectBase.C create mode 100644 src/userobjects/FunctionDPAUserObject.C create mode 100644 src/userobjects/ParkinCoulterDPAUserObject.C rename src/userobjects/{ParkinCoulterBase.C => ParkinCoulterInterface.C} (78%) diff --git a/include/postprocessors/DPAPostprocessor.h b/include/postprocessors/DPAPostprocessor.h new file mode 100644 index 00000000..7ff604d7 --- /dev/null +++ b/include/postprocessors/DPAPostprocessor.h @@ -0,0 +1,34 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#ifdef GSL_ENABLED + +#pragma once + +#include "GeneralPostprocessor.h" + +// forward declarations +class DPAPostprocessor; +class DPAUserObjectBase; + +template <> +InputParameters validParams(); + +class DPAPostprocessor : public GeneralPostprocessor +{ +public: + DPAPostprocessor(const InputParameters & parameters); + virtual void execute() override {} + virtual void initialize() override {} + virtual Real getValue() override; + +protected: + const DPAUserObjectBase & _damage_object; +}; + +#endif diff --git a/include/userobjects/DPAUserObjectBase.h b/include/userobjects/DPAUserObjectBase.h new file mode 100644 index 00000000..7eb452b1 --- /dev/null +++ b/include/userobjects/DPAUserObjectBase.h @@ -0,0 +1,101 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ +#ifdef GSL_ENABLED + +#pragma once + +#include "GeneralUserObject.h" + +class DPAUserObjectBase; + +template <> +InputParameters validParams(); + +class DPAUserObjectBase : public GeneralUserObject +{ +public: + DPAUserObjectBase(const InputParameters & parameters); + + ///@{ get the + Real getDPA() const { return _dpa; } + Real getPartialDPA(unsigned int Z, Real A) const; + ///@} + + void initialize() override {} + +protected: + void prepare(); + bool changed() const; + std::vector getAtomicNumbers() const; + std::vector getMassNumbers() const; + std::vector getNumberFractions() const; + Real getMaxEnergy() const; + + /// accumulated dpa + void accumulateDamage(); + + /// a helper function that sets _ns and checks consistency of A, Z, N + void initAZNHelper(); + + /// a helper function that computes the neutron damage efficiency + Real + neutronDamageEfficiency(unsigned int i, unsigned int j, unsigned int g, unsigned int x) const; + + /// a virtual function computing the integral damage function + virtual Real integralDamageFunction(Real T, unsigned int i, unsigned int j) const = 0; + + /// callback that is executed when composition changed and damage functions must be recomputed + virtual void onCompositionChanged() = 0; + + /// a helper that assigns a unique string to a Z, A pair + std::string zaidHelper(unsigned int Z, Real A) const; + + /// tolerance for recomputing the displacement function + Real _tol = 1e-10; + + /// the computed dose rates + Real _dpa = 0; + + /// the computed dose rate by species; this is a map because presence of (Z,A) can change dynamically + std::map _partial_dpa; + + /// is damage accumulated during a transient or computed for steady state + bool _is_transient_irradiation; + + /// irradiation_time used when dpa is estimated from steady-state calculations + Real _irradiation_time; + + /// the neutron reaction types considered for computing dpa + MultiMooseEnum _neutron_reaction_types; + + /// number of reaction types creating radiation damage + unsigned int _nr; + + ///@{ data used for computing dpa value + std::vector _atomic_numbers; + std::vector _mass_numbers; + std::vector _number_densities; + std::vector _energy_group_boundaries; + std::vector _scalar_flux; + std::vector>> _cross_sections; + ///@} + + /// Q values for each reaction type and isotope + std::vector> _q_values; + + ///@{ the "old" versions of the data; used for determining if disp function update is required + std::vector _atomic_numbers_old; + std::vector _mass_numbers_old; + std::vector _number_densities_old; + ///@} + + /// number of neutron energy groups + unsigned int _ng; +}; + +#endif // GSL_ENABLED diff --git a/include/userobjects/FunctionDPAUserObject.h b/include/userobjects/FunctionDPAUserObject.h new file mode 100644 index 00000000..92af7aa1 --- /dev/null +++ b/include/userobjects/FunctionDPAUserObject.h @@ -0,0 +1,44 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#ifdef GSL_ENABLED + +#pragma once + +#include "DPAUserObjectBase.h" +#include "LinearInterpolation.h" + +class FunctionDPAUserObject; + +template <> +InputParameters validParams(); + +class FunctionDPAUserObject : public DPAUserObjectBase +{ +public: + FunctionDPAUserObject(const InputParameters & parameters); + void finalize() override; + void execute() override; + void initialSetup() override; + void initialize() override {} + +protected: + virtual Real integralDamageFunction(Real T, unsigned int i, unsigned int j) const override; + virtual void onCompositionChanged() override; + + /// the maximum energy step size used for interpolation and integration of integral damage function + Real _max_delta_E; + + /// the damage functions are provided by MOOSE functions + std::vector> _damage_functions; + + /// stores the integral damage functions computed from input as LinearInterpolation objects + std::vector> _integral_damage_functions; +}; + +#endif diff --git a/include/userobjects/ParkinCoulterDPAUserObject.h b/include/userobjects/ParkinCoulterDPAUserObject.h new file mode 100644 index 00000000..ad3a80d3 --- /dev/null +++ b/include/userobjects/ParkinCoulterDPAUserObject.h @@ -0,0 +1,38 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ +#ifdef GSL_ENABLED + +#pragma once + +#include "DPAUserObjectBase.h" +#include "ParkinCoulterInterface.h" + +class ParkinCoulterDPAUserObject; + +template <> +InputParameters validParams(); + +class ParkinCoulterDPAUserObject : public DPAUserObjectBase, public ParkinCoulterInterface +{ +public: + ParkinCoulterDPAUserObject(const InputParameters & parameters); + void finalize() override; + void execute() override; + void initialSetup() override; + +protected: + virtual void initDamageFunctions() override; + virtual std::vector atomicNumbers() const override; + virtual std::vector massNumbers() const override; + virtual std::vector numberFractions() const override; + virtual Real maxEnergy() const override { return getMaxEnergy(); } + virtual Real integralDamageFunction(Real T, unsigned int i, unsigned int j) const override; + virtual void onCompositionChanged() override; +}; + +#endif // GSL_ENABLED diff --git a/include/userobjects/ParkinCoulterBase.h b/include/userobjects/ParkinCoulterInterface.h similarity index 81% rename from include/userobjects/ParkinCoulterBase.h rename to include/userobjects/ParkinCoulterInterface.h index 2889e4d4..7b6cc25d 100644 --- a/include/userobjects/ParkinCoulterBase.h +++ b/include/userobjects/ParkinCoulterInterface.h @@ -16,19 +16,19 @@ // mytrim includes #include -class ParkinCoulterBase; +class ParkinCoulterInterface; class PolyatomicDisplacementFunction; class PolyatomicDamageEnergyFunction; class PolyatomicDisplacementDerivativeFunction; template <> -InputParameters validParams(); +InputParameters validParams(); -class ParkinCoulterBase : public GeneralUserObject +class ParkinCoulterInterface { public: - ParkinCoulterBase(const InputParameters & parameters); - void initialize() override {} + ParkinCoulterInterface(const MooseObject * moose_object); + static InputParameters validParams(); protected: /// recomputes the Polyatomic damage functions @@ -50,6 +50,12 @@ class ParkinCoulterBase : public GeneralUserObject /// computes the polymat object that is used to create PolyatomicDisplacementFunctions std::vector polyMat() const; + /// the MooseObject that inherits from this interface + const MooseObject * _moose_obj; + + /// convenience reference to InputParameters from MooseObject + const InputParameters & _pars; + std::vector> _Ecap; std::unique_ptr _padf; diff --git a/include/userobjects/PolyatomicRecoil.h b/include/userobjects/PolyatomicRecoil.h index a24d22c3..f7651249 100644 --- a/include/userobjects/PolyatomicRecoil.h +++ b/include/userobjects/PolyatomicRecoil.h @@ -9,11 +9,10 @@ #pragma once -#include "ParkinCoulterBase.h" +#include "GeneralUserObject.h" +#include "ParkinCoulterInterface.h" -class PolyatomicRecoil; - -class PolyatomicRecoil : public ParkinCoulterBase +class PolyatomicRecoil : public GeneralUserObject, public ParkinCoulterInterface { public: static InputParameters validParams(); @@ -21,6 +20,7 @@ class PolyatomicRecoil : public ParkinCoulterBase PolyatomicRecoil(const InputParameters & parameters); void finalize() override; void execute() override; + void initialize() override {} protected: virtual void initDamageFunctions() override; diff --git a/src/postprocessors/DPAPostprocessor.C b/src/postprocessors/DPAPostprocessor.C new file mode 100644 index 00000000..570c51eb --- /dev/null +++ b/src/postprocessors/DPAPostprocessor.C @@ -0,0 +1,37 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#ifdef GSL_ENABLED + +#include "DPAPostprocessor.h" +#include "DPAUserObjectBase.h" + +registerMooseObject("MagpieApp", DPAPostprocessor); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addRequiredParam("dpa_object", "The neutronics damage object."); + params.addClassDescription("Retrieves the value of the dpa from a DPAUserObjectBase."); + return params; +} + +DPAPostprocessor::DPAPostprocessor(const InputParameters & params) + : GeneralPostprocessor(params), _damage_object(getUserObject("dpa_object")) +{ +} + +Real +DPAPostprocessor::getValue() +{ + return _damage_object.getDPA(); +} + +#endif diff --git a/src/userobjects/DPAUserObjectBase.C b/src/userobjects/DPAUserObjectBase.C new file mode 100644 index 00000000..4f940b02 --- /dev/null +++ b/src/userobjects/DPAUserObjectBase.C @@ -0,0 +1,406 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#ifdef GSL_ENABLED + +#include "DPAUserObjectBase.h" +#include "PolyatomicDisplacementFunction.h" + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addParam("irradiation_time", "Irradiation time used "); + MultiMooseEnum damage_reaction_types("elastic inelastic"); + params.addParam("damage_reaction_types", + damage_reaction_types, + "The neutron reaction causing radiation damage"); + params.addParam>("Z", "The atomic numbers of all present isotopes"); + params.addParam>("A", "The mass numbers of all present isotopes"); + params.addParam>("number_densities", + "The number densities of all present isotopes"); + params.addParam>("energy_group_boundaries", + "The neutron flux energy group boundaries in units of eV " + "starting with the highest energy group"); + params.addParam>("scalar_flux", + "The values of the neutron scalar flux by energy group " + "starting with the highest energy group"); + params.addParam>>( + "cross_section", + "The values of the cross sections. Each semicolon separated vector contains cross sections " + "for a particular nuclide and reaction type. Each entry must be number of energy groups " + "entries long. One vector must be provided for each " + "nuclide/reaction type combination. The ordering is as follows: if there are reaction types" + "a and b, and nuclides i, k, and l, the ordering will be xs_ai; xs_ak; xs_al; xs_bi; xs_bk, " + "xs_bl"); + params.addParam>>( + "Q", "The Q values by reaction type and then isotope. Assumed zero if not provided."); + params.set("execute_on") = EXEC_TIMESTEP_BEGIN; + params.suppressParameter("execute_on"); + return params; +} + +DPAUserObjectBase::DPAUserObjectBase(const InputParameters & parameters) + : GeneralUserObject(parameters), + _is_transient_irradiation(!isParamValid("irradiation_time")), + _irradiation_time(_is_transient_irradiation ? 0 : getParam("irradiation_time")), + _neutron_reaction_types(getParam("damage_reaction_types")), + _nr(_neutron_reaction_types.size()) +{ + if (_nr == 0) + paramError("damage_reaction_types", "At least one damage mechanism must be provided"); + + // get parameters if they are provided + if (isParamValid("Z")) + _atomic_numbers = getParam>("Z"); + if (isParamValid("A")) + _mass_numbers = getParam>("A"); + if (isParamValid("number_densities")) + _number_densities = getParam>("number_densities"); + if (isParamValid("energy_group_boundaries")) + _energy_group_boundaries = getParam>("energy_group_boundaries"); + if (isParamValid("scalar_flux")) + _scalar_flux = getParam>("scalar_flux"); + + if (isParamValid("cross_section")) + { + std::vector> xs = getParam>>("cross_section"); + + // we know _nr so the length of cross sections should be divisible by _nr + unsigned int n = xs.size(); + if (n % _nr != 0) + paramError("cross_section", + "The number of entries in the cross_section param must be divisible by number of " + "reaction types ", + _nr); + + unsigned int n_nuc = n / _nr; + _cross_sections.resize(_nr); + unsigned int p = 0; + for (unsigned int r = 0; r < _nr; ++r) + { + _cross_sections[r].resize(n_nuc); + for (unsigned int j = 0; j < n_nuc; ++j) + { + _cross_sections[r][j] = xs[p]; + ++p; + } + } + } +} + +void +DPAUserObjectBase::prepare() +{ + // This object support dynamic changes of number densities + // appearance and disappearance of isotopes, and change in + // cross sections. Therefore, this function needs to check + // and redo a lot of things that would usually be done in + // the constructor. + + // checks on scalar flux & energy groups + _ng = _scalar_flux.size(); + if (_ng == 0) + mooseError("The number of provided scalar fluxes is zero. Provide using parameter scalar_flux " + "or set programmatically"); + + if (_energy_group_boundaries.size() != _ng + 1) + mooseError("The size of the energy group boundary is ", + _energy_group_boundaries.size(), + " but it must have exactly number of energy groups plus one entry ", + _ng + 1); + + for (unsigned int j = 1; j < _energy_group_boundaries.size(); ++j) + if (_energy_group_boundaries[j - 1] <= _energy_group_boundaries[j]) + mooseError("Entries of energy_group_boundaries must be strictly decreasing"); + + // check on A, Z, number densities + initAZNHelper(); + + // checks on cross sections + if (_cross_sections.size() != _nr) + mooseError("Leading dimension of _cross_sections is ", + _cross_sections.size(), + " but should be equal to number of damage reaction types ", + _nr); + for (unsigned int i = 0; i < _nr; ++i) + { + if (_cross_sections[i].size() != _atomic_numbers.size()) + mooseError("Second dimension of _cross_sections for index ", + i, + " does not match number of isotope ", + _atomic_numbers.size()); + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + if (_cross_sections[i][j].size() != _ng) + mooseError("Third dimension of _cross_sections for indices ", + i, + " ", + j, + " has ", + _cross_sections[i][j].size(), + " entries when number of groups is ", + _ng); + } + + // get Q values, need to do this here to be sure that number of isotopes is set + if (isParamValid("Q")) + { + _q_values = getParam>>("Q"); + if (_q_values.size() != _nr) + paramError("Q", "Leading dimension must be of the same size as reaction types."); + for (unsigned int j = 0; j < _nr; ++j) + if (_q_values[j].size() != _atomic_numbers.size()) + paramError("Q", "Second dimension must be of same size as isotopes."); + } + else + { + _q_values.resize(_nr); + for (unsigned int j = 0; j < _nr; ++j) + _q_values[j].resize(_atomic_numbers.size()); + } +} + +void +DPAUserObjectBase::initAZNHelper() +{ + unsigned int ns = getAtomicNumbers().size(); + if (ns == 0) + mooseError("The size of the atomic number array is zero. Set Z parameter or set atomic number " + "programmatically"); + + if (_mass_numbers.size() != ns) + mooseError("Size of mass number array does not match size of atomic number array. Size of mass " + "numbers ", + _mass_numbers.size(), + ", size of atomic numbers ", + ns); + + if (_number_densities.size() != ns) + mooseError("Size of mass number array does not match size of atomic number array. Size of mass " + "numbers ", + _number_densities.size(), + ", size of atomic numbers ", + ns); +} + +std::vector +DPAUserObjectBase::getAtomicNumbers() const +{ + // in this kind we need to work a bit since + // VPPs are Reals and we need to return unsigned int + std::vector Zs(_atomic_numbers.size()); + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + { + unsigned int i = static_cast(_atomic_numbers[j]); + if (std::abs(_atomic_numbers[j] - i) > 1e-12) + mooseError("Entry ", j, ":", _atomic_numbers[j], " is not a non-negative integer."); + Zs[j] = i; + } + return Zs; +} + +std::vector +DPAUserObjectBase::getMassNumbers() const +{ + return _mass_numbers; +} + +std::vector +DPAUserObjectBase::getNumberFractions() const +{ + // need to normalize here + Real sum = 0; + for (unsigned int j = 0; j < _number_densities.size(); ++j) + sum += _number_densities[j]; + std::vector nf(_number_densities.size()); + for (unsigned int j = 0; j < _number_densities.size(); ++j) + nf[j] = _number_densities[j] / sum; + return nf; +} + +Real +DPAUserObjectBase::getMaxEnergy() const +{ + // between elastic & inelastic scattering, elastic scattering + // allows the largest recoil energy = gamma * E_neutron + Real max_neutron_energy = _energy_group_boundaries[0]; + + // find maximum gamma + Real min_mass = std::numeric_limits::max(); + for (unsigned int j = 0; j < _mass_numbers.size(); ++j) + if (min_mass > _mass_numbers[j]) + min_mass = _mass_numbers[j]; + Real gamma = 4 * min_mass / (min_mass + 1) / (min_mass + 1); + return gamma * max_neutron_energy; +} + +bool +DPAUserObjectBase::changed() const +{ + // the size could have changed + if (_atomic_numbers.size() != _atomic_numbers_old.size() || + _mass_numbers.size() != _mass_numbers_old.size() || + _number_densities.size() != _number_densities_old.size()) + return true; + + // the size is still the same + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + if (!MooseUtils::absoluteFuzzyEqual(_atomic_numbers[j], _atomic_numbers_old[j], _tol) || + !MooseUtils::absoluteFuzzyEqual(_mass_numbers[j], _mass_numbers_old[j], _tol) || + !MooseUtils::absoluteFuzzyEqual(_number_densities[j], _number_densities_old[j], _tol)) + return true; + + return false; +} + +void +DPAUserObjectBase::accumulateDamage() +{ + if (changed()) + { + initAZNHelper(); + onCompositionChanged(); + + // copy over A, Z, N to A, Z, N_old + _atomic_numbers_old.resize(_atomic_numbers.size()); + _mass_numbers_old.resize(_atomic_numbers.size()); + _number_densities_old.resize(_atomic_numbers.size()); + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + { + _atomic_numbers_old[j] = _atomic_numbers[j]; + _mass_numbers_old[j] = _mass_numbers[j]; + _number_densities_old[j] = _number_densities[j]; + } + } + + // compute total number density + Real total_number_density = 0; + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + total_number_density += _number_densities[j]; + + // damage function have been computed + // so now we can start computing dpa + if (!_is_transient_irradiation) + { + _dpa = 0; + _partial_dpa.clear(); + } + + Real del_t = _is_transient_irradiation ? _dt : _irradiation_time; + + for (unsigned int x = 0; x < _nr; ++x) + { + for (unsigned int g = 0; g < _ng; ++g) + { + Real del_E = _energy_group_boundaries[g] - _energy_group_boundaries[g + 1]; + + for (unsigned int i = 0; i < _atomic_numbers.size(); ++i) + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + { + Real v = del_t / del_E * _scalar_flux[g] * _number_densities[i] / total_number_density * + _cross_sections[x][i][g] * neutronDamageEfficiency(i, j, g, x); + + // total dpa is incremented regardless + _dpa += v; + + // partial dpa: if the entry exists, add to it, otherwise create it + const auto & p = _partial_dpa.find(zaidHelper(_atomic_numbers[i], _mass_numbers[i])); + if (p == _partial_dpa.end()) + _partial_dpa[zaidHelper(_atomic_numbers[i], _mass_numbers[i])] = v; + else + p->second += v; + } + } + } +} + +Real +DPAUserObjectBase::neutronDamageEfficiency(unsigned int i, + unsigned int j, + unsigned int g, + unsigned int x) const +{ + std::vector points; + std::vector weights; + PolyatomicDisplacementFunction::gslQuadRule(100, points, weights); + + if (_neutron_reaction_types[x] == "elastic") + { + Real gamma = 4 * _mass_numbers[i] / (_mass_numbers[i] + 1) / (_mass_numbers[i] + 1); + Real from_E = _energy_group_boundaries[g + 1]; + Real to_E = _energy_group_boundaries[g]; + Real delta_E = to_E - from_E; + + // integral over E + Real integral = 0; + for (unsigned int qp = 0; qp < points.size(); ++qp) + { + Real energy = 0.5 * (points[qp] + 1) * delta_E + from_E; + Real w = 0.5 * weights[qp] * delta_E; + integral += w * integralDamageFunction(gamma * energy, i, j) / energy; + } + + return integral / gamma; + } + else if (_neutron_reaction_types[x] == "inelastic") + { + // warn if q value does not make sense + if (_q_values[x][i] >= 0) + mooseDoOnce(mooseWarning("Q value is greater or equal to zero for inelastic scattering.")); + + Real d = std::abs(_q_values[x][i]) * (_mass_numbers[i] + 1) / _mass_numbers[i]; + Real gamma = 4 * _mass_numbers[i] / (_mass_numbers[i] + 1) / (_mass_numbers[i] + 1); + Real from_E = _energy_group_boundaries[g + 1]; + Real to_E = _energy_group_boundaries[g]; + Real delta_E = to_E - from_E; + + // integral over E + Real integral = 0; + for (unsigned int qp = 0; qp < points.size(); ++qp) + { + Real energy = 0.5 * (points[qp] + 1) * delta_E + from_E; + Real Delta = d / energy; + + // make sure threshold energy is honored + if (Delta > 1) + continue; + + Real sr = std::sqrt(1 - Delta); + Real Tmin = gamma / 2 * energy * (1 - sr - 0.5 * Delta); + Real Tmax = gamma / 2 * energy * (1 + sr - 0.5 * Delta); + Real w = 0.5 * weights[qp] * delta_E; + integral += w * (integralDamageFunction(Tmax, i, j) - integralDamageFunction(Tmin, i, j)) / + (energy * sr); + } + + return integral / gamma; + } + + mooseError("Neutron reaction type not recognized. Should never get here."); + return 0; +} + +std::string +DPAUserObjectBase::zaidHelper(unsigned int Z, Real A) const +{ + std::stringstream ss; + ss << std::setprecision(4) << Z << A; + return ss.str(); +} + +Real +DPAUserObjectBase::getPartialDPA(unsigned int Z, Real A) const +{ + const auto & p = _partial_dpa.find(zaidHelper(Z, A)); + if (p == _partial_dpa.end()) + return 0; + return p->second; +} + +#endif diff --git a/src/userobjects/FunctionDPAUserObject.C b/src/userobjects/FunctionDPAUserObject.C new file mode 100644 index 00000000..a56e66d9 --- /dev/null +++ b/src/userobjects/FunctionDPAUserObject.C @@ -0,0 +1,151 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#ifdef GSL_ENABLED + +#include "FunctionDPAUserObject.h" +#include "PolyatomicDisplacementFunction.h" + +registerMooseObject("MagpieApp", FunctionDPAUserObject); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addRequiredParam>>( + "damage_functions", "Damage functions for each combinations of projectiles and targets."); + params.addParam( + "max_energy_step_size", + 100, + "The maximum energy step size used for integration and interpolation. Default is 100 eV."); + params.addClassDescription("Computes the dose in dpa from composition, cross section, damage " + "type, and neutron flux. The damage functions are provided as MOOSE " + "functions where energy is provided via the time arg slot."); + return params; +} + +FunctionDPAUserObject::FunctionDPAUserObject(const InputParameters & parameters) + : DPAUserObjectBase(parameters), _max_delta_E(getParam("max_energy_step_size")) +{ + // get the damage functions + std::vector> nm = + getParam>>("damage_functions"); + _damage_functions.resize(nm.size()); + for (unsigned int j = 0; j < nm.size(); ++j) + { + _damage_functions[j].resize(nm[j].size()); + for (unsigned int i = 0; i < nm[j].size(); ++i) + _damage_functions[j][i] = &getFunctionByName(nm[j][i]); + } +} + +void +FunctionDPAUserObject::initialSetup() +{ + prepare(); + + // make sure the dimension of _damage_functions is correct + if (_damage_functions.size() != _atomic_numbers.size()) + paramError("damage_functions", + "Size of leading dimension ", + _damage_functions.size(), + " is not identical to number of isotopes ", + _atomic_numbers.size()); + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + if (_damage_functions[j].size() != _atomic_numbers.size()) + paramError("damage_functions", + "Size of entry ", + j, + " is ", + _damage_functions[j].size(), + " which is not identical to number of isotopes ", + _atomic_numbers.size()); +} + +void +FunctionDPAUserObject::execute() +{ + accumulateDamage(); +} + +void +FunctionDPAUserObject::onCompositionChanged() +{ + // compute the integral damage functions + _integral_damage_functions.resize(_atomic_numbers.size()); + for (unsigned int i = 0; i < _atomic_numbers.size(); ++i) + { + // resize integral damage function + _integral_damage_functions[i].resize(_atomic_numbers.size()); + + for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) + { + // loop over all energies + Real energy = 0; + Real max_energy = getMaxEnergy(); + std::vector energies(1); + std::vector integral(1); + // this is the energy time step, initially it must be + // small but then can up to _max_delta_E + Real deltaE = 1e-5; + bool keep_going = true; + while (keep_going) + { + Real energy_old = energy; + energy += deltaE; + if (energy >= max_energy) + { + energy = max_energy; + keep_going = false; + } + + // incremental integration from energy_old to energy + std::vector points; + std::vector weights; + PolyatomicDisplacementFunction::gslQuadRule(100, points, weights); + + Real tint = 0; + for (unsigned int qp = 0; qp < points.size(); ++qp) + { + Real e = 0.5 * (points[qp] + 1) * deltaE + energy_old; + Real w = 0.5 * weights[qp] * deltaE; + tint += w * _damage_functions[i][j]->value(e, Point()); + } + + // update energies and integral vectors + energies.push_back(energy); + Real old_integral = integral.back(); + integral.push_back(tint + old_integral); + + // now do an adjustment of the timestep but only up to _max_delta_E + // the factor of 1.01 was found to be well converged for verification_2.i + Real proposed_deltaE = deltaE * 1.01; + if (proposed_deltaE < _max_delta_E) + deltaE = proposed_deltaE; + } + + // now the linear interpolation object is constructed and stored in integral damage function + // var + _integral_damage_functions[i][j] = LinearInterpolation(energies, integral); + } + } +} + +Real +FunctionDPAUserObject::integralDamageFunction(Real T, unsigned int i, unsigned int j) const +{ + return _integral_damage_functions[i][j].sample(T); +} + +void +FunctionDPAUserObject::finalize() +{ +} + +#endif diff --git a/src/userobjects/ParkinCoulterDPAUserObject.C b/src/userobjects/ParkinCoulterDPAUserObject.C new file mode 100644 index 00000000..7bbbbb4e --- /dev/null +++ b/src/userobjects/ParkinCoulterDPAUserObject.C @@ -0,0 +1,94 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ +#ifdef GSL_ENABLED + +#include "ParkinCoulterDPAUserObject.h" +#include "PolyatomicRecoil.h" +#include "PolyatomicDisplacementFunction.h" +#include "PolyatomicDamageEnergyFunction.h" +#include "PolyatomicDisplacementDerivativeFunction.h" +#include "MooseMesh.h" + +// various other includes +#include +#include + +registerMooseObject("MagpieApp", ParkinCoulterDPAUserObject); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params += ParkinCoulterInterface::validParams(); + params.addClassDescription( + "Computes the dose in dpa from composition, cross section, damage type, and neutron flux for " + "polyatomic materials using Parkin-Coulter's method."); + return params; +} + +ParkinCoulterDPAUserObject::ParkinCoulterDPAUserObject(const InputParameters & parameters) + : DPAUserObjectBase(parameters), ParkinCoulterInterface(this) +{ +} + +void +ParkinCoulterDPAUserObject::initialSetup() +{ + prepare(); +} + +std::vector +ParkinCoulterDPAUserObject::atomicNumbers() const +{ + return getAtomicNumbers(); +} + +std::vector +ParkinCoulterDPAUserObject::massNumbers() const +{ + return getMassNumbers(); +} + +std::vector +ParkinCoulterDPAUserObject::numberFractions() const +{ + return getNumberFractions(); +} + +void +ParkinCoulterDPAUserObject::initDamageFunctions() +{ + _padf = libmesh_make_unique(polyMat(), NET, _Ecap); +} + +void +ParkinCoulterDPAUserObject::execute() +{ + accumulateDamage(); +} + +void +ParkinCoulterDPAUserObject::onCompositionChanged() +{ + computeDamageFunctions(); + _padf->computeDisplacementFunctionIntegral(); +} + +Real +ParkinCoulterDPAUserObject::integralDamageFunction(Real T, unsigned int i, unsigned int j) const +{ + return _padf->linearInterpolationIntegralDamageFunction(T, i, j); +} + +void +ParkinCoulterDPAUserObject::finalize() +{ +} + +#endif diff --git a/src/userobjects/ParkinCoulterBase.C b/src/userobjects/ParkinCoulterInterface.C similarity index 78% rename from src/userobjects/ParkinCoulterBase.C rename to src/userobjects/ParkinCoulterInterface.C index 550cae9b..62743399 100644 --- a/src/userobjects/ParkinCoulterBase.C +++ b/src/userobjects/ParkinCoulterInterface.C @@ -7,7 +7,7 @@ /**********************************************************************/ #ifdef GSL_ENABLED -#include "ParkinCoulterBase.h" +#include "ParkinCoulterInterface.h" #include "PolyatomicDisplacementFunction.h" #include "PolyatomicDamageEnergyFunction.h" #include "PolyatomicDisplacementDerivativeFunction.h" @@ -16,11 +16,10 @@ // mytrim includes #include -template <> InputParameters -validParams() +ParkinCoulterInterface::validParams() { - InputParameters params = validParams(); + InputParameters params = emptyInputParameters(); params.addRequiredParam>("displacement_thresholds", "Dispacement thresholds"); params.addParam>("lattice_binding_energies", "Lattice binding energies"); params.addParam>>( @@ -37,32 +36,28 @@ validParams() "logarithmic_energy_spacing", "logarithmic_energy_spacing > 1", "Spacing of the energy points En in log space energy_spacing = E_{n+1} / En"); - - // this should only be run once - params.set("execute_on") = EXEC_TIMESTEP_BEGIN; - params.suppressParameter("execute_on"); return params; } -ParkinCoulterBase::ParkinCoulterBase(const InputParameters & parameters) - : GeneralUserObject(parameters) +ParkinCoulterInterface::ParkinCoulterInterface(const MooseObject * moose_object) + : _moose_obj(moose_object), _pars(moose_object->parameters()) { _Ecap = {{}}; - if (isParamValid("Ecap")) - _Ecap = getParam>>("Ecap"); + if (_pars.isParamValid("Ecap")) + _Ecap = _pars.get>>("Ecap"); } std::vector -ParkinCoulterBase::polyMat() const +ParkinCoulterInterface::polyMat() const { std::vector poly_mat; std::vector atomic_numbers = atomicNumbers(); std::vector mass_numbers = massNumbers(); std::vector N = numberFractions(); - std::vector threshold = getParam>("displacement_thresholds"); + std::vector threshold = _pars.get>("displacement_thresholds"); std::vector bind; - if (isParamValid("lattice_binding_energies")) - bind = getParam>("lattice_binding_energies"); + if (_pars.isParamValid("lattice_binding_energies")) + bind = _pars.get>("lattice_binding_energies"); else bind.assign(atomic_numbers.size(), 0.0); @@ -87,16 +82,16 @@ ParkinCoulterBase::polyMat() const } void -ParkinCoulterBase::computeDamageFunctions() +ParkinCoulterInterface::computeDamageFunctions() { // callback for allocating damage functions initDamageFunctions(); Real energy = _padf->minEnergy(); Real Emax = maxEnergy(); - Real threshold = getParam("uniform_energy_spacing_threshold"); - Real dE = getParam("uniform_energy_spacing"); - Real logdE = getParam("logarithmic_energy_spacing"); + Real threshold = _pars.get("uniform_energy_spacing_threshold"); + Real dE = _pars.get("uniform_energy_spacing"); + Real logdE = _pars.get("logarithmic_energy_spacing"); for (;;) // while (energy <= Emax) { diff --git a/src/userobjects/PolyatomicRecoil.C b/src/userobjects/PolyatomicRecoil.C index c74ad19d..ecf97882 100644 --- a/src/userobjects/PolyatomicRecoil.C +++ b/src/userobjects/PolyatomicRecoil.C @@ -21,7 +21,8 @@ registerMooseObject("MagpieApp", PolyatomicRecoil); InputParameters PolyatomicRecoil::validParams() { - InputParameters params = ParkinCoulterBase::validParams(); + InputParameters params = GeneralUserObject::validParams(); + params += ParkinCoulterInterface::validParams(); params.addRequiredParam>("Z", "Atomic numbers"); params.addRequiredParam>("A", "Mass numbers"); params.addRequiredParam>("number_fraction", "Number fractions"); @@ -48,7 +49,7 @@ PolyatomicRecoil::validParams() } PolyatomicRecoil::PolyatomicRecoil(const InputParameters & parameters) - : ParkinCoulterBase(parameters) + : GeneralUserObject(parameters), ParkinCoulterInterface(this) { }