diff --git a/include/userobjects/ParkinCoulterBase.h b/include/userobjects/ParkinCoulterBase.h new file mode 100644 index 00000000..2889e4d4 --- /dev/null +++ b/include/userobjects/ParkinCoulterBase.h @@ -0,0 +1,59 @@ +/**********************************************************************/ +/* 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" + +#include "PolyatomicDisplacementFunctionBase.h" + +// mytrim includes +#include + +class ParkinCoulterBase; +class PolyatomicDisplacementFunction; +class PolyatomicDamageEnergyFunction; +class PolyatomicDisplacementDerivativeFunction; + +template <> +InputParameters validParams(); + +class ParkinCoulterBase : public GeneralUserObject +{ +public: + ParkinCoulterBase(const InputParameters & parameters); + void initialize() override {} + +protected: + /// recomputes the Polyatomic damage functions + void computeDamageFunctions(); + + /// this function is called in computeDamageFunctions and is intended to allow + /// allocation of _padf &| _padf_derivative + virtual void initDamageFunctions() = 0; + + ///@{ these functions provide Z, A, and number fraction and must be overidden + virtual std::vector atomicNumbers() const = 0; + virtual std::vector massNumbers() const = 0; + virtual std::vector numberFractions() const = 0; + ///@} + + /// this function provides the maximum energy to integrate PC eqs to + virtual Real maxEnergy() const = 0; + + /// computes the polymat object that is used to create PolyatomicDisplacementFunctions + std::vector polyMat() const; + + std::vector> _Ecap; + + std::unique_ptr _padf; + std::unique_ptr _padf_derivative; +}; + +#endif // GSL_ENABLED diff --git a/include/userobjects/PolyatomicRecoil.h b/include/userobjects/PolyatomicRecoil.h index 555caf70..d6664cab 100644 --- a/include/userobjects/PolyatomicRecoil.h +++ b/include/userobjects/PolyatomicRecoil.h @@ -9,31 +9,26 @@ #pragma once -#include "GeneralUserObject.h" +#include "ParkinCoulterBase.h" class PolyatomicRecoil; -class PolyatomicDisplacementFunction; -class PolyatomicDamageEnergyFunction; -class PolyatomicDisplacementDerivativeFunction; template <> InputParameters validParams(); -class PolyatomicRecoil : public GeneralUserObject +class PolyatomicRecoil : public ParkinCoulterBase { public: PolyatomicRecoil(const InputParameters & parameters); - - void execute() override; - void initialize() override {} void finalize() override; + void execute() override; protected: - std::vector _atomic_numbers; - std::vector _mass_numbers; - - std::unique_ptr _padf; - std::unique_ptr _padf_derivative; + 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 getParam("Emax"); } }; #endif // GSL_ENABLED diff --git a/src/userobjects/ParkinCoulterBase.C b/src/userobjects/ParkinCoulterBase.C new file mode 100644 index 00000000..550cae9b --- /dev/null +++ b/src/userobjects/ParkinCoulterBase.C @@ -0,0 +1,119 @@ +/**********************************************************************/ +/* 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 "ParkinCoulterBase.h" +#include "PolyatomicDisplacementFunction.h" +#include "PolyatomicDamageEnergyFunction.h" +#include "PolyatomicDisplacementDerivativeFunction.h" +#include "MooseMesh.h" + +// mytrim includes +#include + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addRequiredParam>("displacement_thresholds", "Dispacement thresholds"); + params.addParam>("lattice_binding_energies", "Lattice binding energies"); + params.addParam>>( + "Ecap", "Capture energy Ecap_ij of species i being trapped in j site"); + params.addRangeCheckedParam("uniform_energy_spacing_threshold", + 10, + "uniform_energy_spacing_threshold >= 0", + "Threshold below which energy points are spaced uniformly."); + params.addRangeCheckedParam("uniform_energy_spacing", + 0.25, + "uniform_energy_spacing > 0", + "Uniform energy spacing below the threshold"); + params.addRequiredRangeCheckedParam( + "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) +{ + _Ecap = {{}}; + if (isParamValid("Ecap")) + _Ecap = getParam>>("Ecap"); +} + +std::vector +ParkinCoulterBase::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 bind; + if (isParamValid("lattice_binding_energies")) + bind = getParam>("lattice_binding_energies"); + else + bind.assign(atomic_numbers.size(), 0.0); + + // perform some checks + if (atomic_numbers.size() != mass_numbers.size() || atomic_numbers.size() != N.size() || + atomic_numbers.size() != threshold.size() || atomic_numbers.size() != bind.size()) + mooseError("Size mismatch for at least one parameter array. Z, A, number_fraction, " + "displacement_thresholds and lattice_binding_energies" + "must all have the same length."); + + for (unsigned int j = 0; j < atomic_numbers.size(); ++j) + { + MyTRIM_NS::Element element; + element._Z = atomic_numbers[j]; + element._m = mass_numbers[j]; + element._t = N[j]; + element._Edisp = threshold[j]; + element._Elbind = bind[j]; + poly_mat.push_back(element); + } + return poly_mat; +} + +void +ParkinCoulterBase::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"); + + for (;;) // while (energy <= Emax) + { + energy = energy < threshold ? energy + dE : energy * logdE; + if (energy > Emax) + { + _padf->advanceDisplacements(Emax); + break; + } + + // increment displacements for value of energy + _padf->advanceDisplacements(energy); + } + + if (_padf_derivative) + for (unsigned int n = 1; n < _padf->nEnergySteps(); ++n) + _padf_derivative->advanceDisplacements(_padf->energyPoint(n)); +} + +#endif diff --git a/src/userobjects/PolyatomicRecoil.C b/src/userobjects/PolyatomicRecoil.C index b9eae0fd..fa9c309c 100644 --- a/src/userobjects/PolyatomicRecoil.C +++ b/src/userobjects/PolyatomicRecoil.C @@ -22,28 +22,10 @@ template <> InputParameters validParams() { - InputParameters params = validParams(); + InputParameters params = validParams(); params.addRequiredParam>("Z", "Atomic numbers"); params.addRequiredParam>("A", "Mass numbers"); params.addRequiredParam>("number_fraction", "Number fractions"); - params.addRequiredParam>("displacement_thresholds", "Dispacement thresholds"); - params.addParam>("lattice_binding_energies", "Lattice binding energies"); - params.addParam>>( - "Ecap", "Capture energy Ecap_ij of species i being trapped in j site"); - params.addRangeCheckedParam("uniform_energy_spacing_threshold", - 10, - "uniform_energy_spacing_threshold >= 0", - "Threshold below which energy points are spaced uniformly."); - params.addRangeCheckedParam("uniform_energy_spacing", - 0.25, - "uniform_energy_spacing > 0", - "Uniform energy spacing below the threshold"); - params.addRequiredRangeCheckedParam( - "logarithmic_energy_spacing", - "logarithmic_energy_spacing > 1", - "Spacing of the energy points En in log space energy_spacing = E_{n+1} / En"); - params.addRequiredRangeCheckedParam( - "Emax", "Emax > 0", "Maximum desired energy to which displacement functions are computed"); MooseEnum nrt_damage_types("TOTAL NET ENERGY NET_DERIVATIVE", "TOTAL"); params.addParam( "damage_type", @@ -55,48 +37,43 @@ validParams() "NET_DERIVATIVE: derivative of NET w.r.t. the partial number fractions."); params.addParam("displacement_file_base", "The output file base for displacement function in csv format."); + params.addRequiredRangeCheckedParam( + "Emax", "Emax > 0", "Maximum desired energy to which displacement functions are computed"); params.addClassDescription( "PolyatomicRecoil allows computation of total and net displacement functions," "damage energy functions, and the derivative of the net displacement functions w.r.t. number " "fractions."); + params.set("execute_on") = EXEC_INITIAL; + params.suppressParameter("execute_on"); return params; } PolyatomicRecoil::PolyatomicRecoil(const InputParameters & parameters) - : GeneralUserObject(parameters), - _atomic_numbers(getParam>("Z")), - _mass_numbers(getParam>("A")) + : ParkinCoulterBase(parameters) { - std::vector N = getParam>("number_fraction"); - std::vector threshold = getParam>("displacement_thresholds"); - std::vector bind; - if (isParamValid("lattice_binding_energies")) - bind = getParam>("lattice_binding_energies"); - else - bind.assign(_atomic_numbers.size(), 0.0); +} - std::vector> Ecap = {{}}; - if (isParamValid("Ecap")) - Ecap = getParam>>("Ecap"); +std::vector +PolyatomicRecoil::atomicNumbers() const +{ + return getParam>("Z"); +} - if (_atomic_numbers.size() != _mass_numbers.size() || _atomic_numbers.size() != N.size() || - _atomic_numbers.size() != threshold.size() || _atomic_numbers.size() != bind.size()) - mooseError("Size mismatch for at least one parameter array. Z, A, number_fraction, " - "displacement_thresholds and lattice_binding_energies" - "must all have the same length."); +std::vector +PolyatomicRecoil::massNumbers() const +{ + return getParam>("A"); +} - std::vector poly_mat; - for (unsigned int j = 0; j < _atomic_numbers.size(); ++j) - { - MyTRIM_NS::Element element; - element._Z = _atomic_numbers[j]; - element._m = _mass_numbers[j]; - element._t = N[j]; - element._Edisp = threshold[j]; - element._Elbind = bind[j]; - poly_mat.push_back(element); - } +std::vector +PolyatomicRecoil::numberFractions() const +{ + return getParam>("number_fraction"); +} +void +PolyatomicRecoil::initDamageFunctions() +{ // set the displacement function type nrt_type type = TOTAL; if (getParam("damage_type") == "NET") @@ -107,50 +84,29 @@ PolyatomicRecoil::PolyatomicRecoil(const InputParameters & parameters) type = NET_DERIVATIVE; if (type == ENERGY) - _padf = libmesh_make_unique(poly_mat, type, Ecap); + _padf = libmesh_make_unique(polyMat(), type, _Ecap); else if (type == NET_DERIVATIVE) { - _padf = libmesh_make_unique(poly_mat, NET, Ecap); + _padf = libmesh_make_unique(polyMat(), NET, _Ecap); _padf_derivative = libmesh_make_unique( - poly_mat, type, dynamic_cast(_padf.get()), Ecap); + polyMat(), type, dynamic_cast(_padf.get()), _Ecap); } else - _padf = libmesh_make_unique(poly_mat, type, Ecap); - - if (getParam("Emax") < getParam("uniform_energy_spacing_threshold")) - mooseError("Emax must be larger than uniform_energy_spacing_threshold."); + _padf = libmesh_make_unique(polyMat(), type, _Ecap); } void PolyatomicRecoil::execute() { - Real energy = _padf->minEnergy(); - Real Emax = getParam("Emax"); - Real threshold = getParam("uniform_energy_spacing_threshold"); - Real dE = getParam("uniform_energy_spacing"); - Real logdE = getParam("logarithmic_energy_spacing"); - - for (;;) // while (energy <= Emax) - { - energy = energy < threshold ? energy + dE : energy * logdE; - if (energy > Emax) - { - _padf->advanceDisplacements(Emax); - break; - } - - // increment displacements for value of energy - _padf->advanceDisplacements(energy); - } - - if (_padf_derivative) - for (unsigned int n = 1; n < _padf->nEnergySteps(); ++n) - _padf_derivative->advanceDisplacements(_padf->energyPoint(n)); + computeDamageFunctions(); } void PolyatomicRecoil::finalize() { + std::vector atomic_numbers = atomicNumbers(); + std::vector mass_numbers = massNumbers(); + // displacement functions if (isParamValid("displacement_file_base")) { @@ -169,10 +125,10 @@ PolyatomicRecoil::finalize() for (unsigned int derivative = 0; derivative < _padf_derivative->nSpecies(); ++derivative) { unsigned int projectile_zaid = - 1000 * _atomic_numbers[projectile] + _mass_numbers[projectile]; - unsigned int target_zaid = 1000 * _atomic_numbers[target] + _mass_numbers[target]; + 1000 * atomic_numbers[projectile] + mass_numbers[projectile]; + unsigned int target_zaid = 1000 * atomic_numbers[target] + mass_numbers[target]; unsigned int derivative_zaid = - 1000 * _atomic_numbers[derivative] + _mass_numbers[derivative]; + 1000 * atomic_numbers[derivative] + mass_numbers[derivative]; displacement_file << ", d(" << projectile_zaid << "->" << target_zaid << ") / d(" << derivative_zaid << ")"; } @@ -183,15 +139,15 @@ PolyatomicRecoil::finalize() for (unsigned int target = 0; target < _padf->nSpecies(); ++target) { unsigned int projectile_zaid = - 1000 * _atomic_numbers[projectile] + _mass_numbers[projectile]; - unsigned int target_zaid = 1000 * _atomic_numbers[target] + _mass_numbers[target]; + 1000 * atomic_numbers[projectile] + mass_numbers[projectile]; + unsigned int target_zaid = 1000 * atomic_numbers[target] + mass_numbers[target]; displacement_file << "," << projectile_zaid << "->" << target_zaid; } else if (energy_function) for (unsigned int projectile = 0; projectile < _padf->nSpecies(); ++projectile) { unsigned int projectile_zaid = - 1000 * _atomic_numbers[projectile] + _mass_numbers[projectile]; + 1000 * atomic_numbers[projectile] + mass_numbers[projectile]; displacement_file << "," << projectile_zaid; } else