Refactor PolyatomicRecoil by adding base class (idaholab#410)

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 <mytrim/element.h>
+
+class ParkinCoulterBase;
+class PolyatomicDisplacementFunction;
+class PolyatomicDamageEnergyFunction;
+class PolyatomicDisplacementDerivativeFunction;
+
+template <>
+InputParameters validParams<ParkinCoulterBase>();
+
+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<unsigned int> atomicNumbers() const = 0;
+  virtual std::vector<Real> massNumbers() const = 0;
+  virtual std::vector<Real> 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<MyTRIM_NS::Element> polyMat() const;
+
+  std::vector<std::vector<Real>> _Ecap;
+
+  std::unique_ptr<PolyatomicDisplacementFunctionBase> _padf;
+  std::unique_ptr<PolyatomicDisplacementDerivativeFunction> _padf_derivative;
+};
+
+#endif // GSL_ENABLED

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<PolyatomicRecoil>();
 
-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<unsigned int> _atomic_numbers;
-  std::vector<Real> _mass_numbers;
-
-  std::unique_ptr<PolyatomicDisplacementFunctionBase> _padf;
-  std::unique_ptr<PolyatomicDisplacementDerivativeFunction> _padf_derivative;
+  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 getParam<Real>("Emax"); }
 };
 
 #endif // GSL_ENABLED

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 <mytrim/element.h>
+
+template <>
+InputParameters
+validParams<ParkinCoulterBase>()
+{
+  InputParameters params = validParams<GeneralUserObject>();
+  params.addRequiredParam<std::vector<Real>>("displacement_thresholds", "Dispacement thresholds");
+  params.addParam<std::vector<Real>>("lattice_binding_energies", "Lattice binding energies");
+  params.addParam<std::vector<std::vector<Real>>>(
+      "Ecap", "Capture energy Ecap_ij of species i being trapped in j site");
+  params.addRangeCheckedParam<Real>("uniform_energy_spacing_threshold",
+                                    10,
+                                    "uniform_energy_spacing_threshold >= 0",
+                                    "Threshold below which energy points are spaced uniformly.");
+  params.addRangeCheckedParam<Real>("uniform_energy_spacing",
+                                    0.25,
+                                    "uniform_energy_spacing > 0",
+                                    "Uniform energy spacing below the threshold");
+  params.addRequiredRangeCheckedParam<Real>(
+      "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<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
+  params.suppressParameter<ExecFlagEnum>("execute_on");
+  return params;
+}
+
+ParkinCoulterBase::ParkinCoulterBase(const InputParameters & parameters)
+  : GeneralUserObject(parameters)
+{
+  _Ecap = {{}};
+  if (isParamValid("Ecap"))
+    _Ecap = getParam<std::vector<std::vector<Real>>>("Ecap");
+}
+
+std::vector<MyTRIM_NS::Element>
+ParkinCoulterBase::polyMat() const
+{
+  std::vector<MyTRIM_NS::Element> poly_mat;
+  std::vector<unsigned int> atomic_numbers = atomicNumbers();
+  std::vector<Real> mass_numbers = massNumbers();
+  std::vector<Real> N = numberFractions();
+  std::vector<Real> threshold = getParam<std::vector<Real>>("displacement_thresholds");
+  std::vector<Real> bind;
+  if (isParamValid("lattice_binding_energies"))
+    bind = getParam<std::vector<Real>>("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<Real>("uniform_energy_spacing_threshold");
+  Real dE = getParam<Real>("uniform_energy_spacing");
+  Real logdE = getParam<Real>("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