Implement dpa userobject; inherit from NeutronDamageInterface to allo…

…w different damage functions to be plugged in (idaholab#410)

include/userobjects/FunctionDPAUserObject.h

@@ -0,0 +1,45 @@
+/**********************************************************************/
+/*                     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 "NeutronDamageInterface.h"
+#include "LinearInterpolation.h"
+
+class FunctionDPAUserObject;
+
+template <>
+InputParameters validParams<FunctionDPAUserObject>();
+
+class FunctionDPAUserObject : public GeneralUserObject, public NeutronDamageInterface
+{
+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<std::vector<const Function *>> _damage_functions;
+
+  /// stores the integral damage functions computed from input as LinearInterpolation objects
+  std::vector<std::vector<LinearInterpolation>> _integral_damage_functions;
+};
+
+#endif

include/userobjects/NeutronDamageInterface.h

@@ -0,0 +1,115 @@
+/**********************************************************************/
+/*                     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 "MooseTypes.h"
+#include "ExecFlagEnum.h"
+#include "MooseEnum.h"
+#include "InputParameters.h"
+
+class NeutronDamageInterface;
+
+template <>
+InputParameters validParams<NeutronDamageInterface>();
+
+class NeutronDamageInterface
+{
+public:
+  NeutronDamageInterface(const MooseObject * moose_object);
+  static InputParameters validParams();
+
+  ///@{ get the
+  Real getDPA() const { return _dpa; }
+  Real getPartialDPA(unsigned int Z, Real A) const;
+  ///@}
+
+protected:
+  void prepare();
+  bool changed() const;
+  std::vector<unsigned int> getAtomicNumbers() const;
+  std::vector<Real> getMassNumbers() const;
+  std::vector<Real> 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<std::string, Real> _partial_dpa;
+
+  /// the MooseObject that inherits from this interface
+  const MooseObject * _moose_obj;
+
+  /// convenience reference to InputParameters from MooseObject
+  const InputParameters & _pars;
+
+  /// Reference to FEProblemBase
+  FEProblemBase & _nxsi_fe_problem;
+
+  /// 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;
+
+  /// convenience copy of this parameter
+  std::vector<std::vector<std::string>> _xs_names;
+
+  ///@{ data used for computing dpa value
+  const VectorPostprocessorValue & _atomic_numbers;
+  const VectorPostprocessorValue & _mass_numbers;
+  const VectorPostprocessorValue & _number_densities;
+  const VectorPostprocessorValue & _energy_group_boundaries;
+  const VectorPostprocessorValue & _scalar_flux;
+  std::vector<std::vector<const VectorPostprocessorValue *>> _cross_sections;
+  ///@}
+
+  /// Q values for each reaction type and isotope
+  std::vector<std::vector<Real>> _q_values;
+
+  ///@{ the "old" versions of the data; used for determining if disp function update is required
+  std::vector<Real> _atomic_numbers_old;
+  std::vector<Real> _mass_numbers_old;
+  std::vector<Real> _number_densities_old;
+  ///@}
+
+  /// number of neutron energy groups
+  unsigned int _ng;
+
+  /// number of reaction types creating radiation damage
+  unsigned int _nr;
+};
+
+#endif // GSL_ENABLED

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 "ParkinCoulterBase.h"
+#include "NeutronDamageInterface.h"
+
+class ParkinCoulterDPAUserObject;
+
+template <>
+InputParameters validParams<ParkinCoulterDPAUserObject>();
+
+class ParkinCoulterDPAUserObject : public ParkinCoulterBase, public NeutronDamageInterface
+{
+public:
+  ParkinCoulterDPAUserObject(const InputParameters & parameters);
+  void finalize() override;
+  void execute() override;
+  void initialSetup() override;
+
+protected:
+  virtual void initDamageFunctions() override;
+  virtual std::vector<unsigned int> atomicNumbers() const override;
+  virtual std::vector<Real> massNumbers() const override;
+  virtual std::vector<Real> 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

src/userobjects/FunctionDPAUserObject.C

@@ -0,0 +1,156 @@
+/**********************************************************************/
+/*                     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<FunctionDPAUserObject>()
+{
+  InputParameters params = validParams<GeneralUserObject>();
+  params += NeutronDamageInterface::validParams();
+  params.addRequiredParam<std::vector<std::vector<FunctionName>>>(
+      "damage_functions", "Damage functions for each combinations of projectiles and targets.");
+  params.addParam<Real>(
+      "max_energy_step_size",
+      100,
+      "The maximum energy step size used for integration and interpolation. Default is 100 eV.");
+  params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_END;
+  params.suppressParameter<ExecFlagEnum>("execute_on");
+  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)
+  : GeneralUserObject(parameters),
+    NeutronDamageInterface(this),
+    _max_delta_E(getParam<Real>("max_energy_step_size"))
+{
+  // get the damage functions
+  std::vector<std::vector<FunctionName>> nm =
+      getParam<std::vector<std::vector<FunctionName>>>("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<Real> energies(1);
+      std::vector<Real> 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<Real> points;
+        std::vector<Real> 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