PKASurfaceFluxGenerator

include/userobjects/PKASurfaceFluxGenerator.h

@@ -0,0 +1,76 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#pragma once
+
+#include "PKAGeneratorBase.h"
+
+/**
+ * Starts PKAs at a fixed point in a fixed direction
+ */
+class PKASurfaceFluxGenerator : public PKAGeneratorBase
+{
+public:
+  static InputParameters validParams();
+
+  PKASurfaceFluxGenerator(const InputParameters & parameters);
+
+  virtual void appendPKAs(std::vector<MyTRIM_NS::IonBase> &,
+                          const MyTRIMRasterizer::PKAParameters &,
+                          const MyTRIMRasterizer::AveragedData &) const;
+
+protected:
+  /// provides a mean to override the angular distribution of the PKAs in derived class
+  virtual void setDirection(MyTRIM_NS::IonBase & ion) const;
+
+  /// the direction along which the PKAs start
+  const RealVectorValue _direction;
+
+  /// Uses point locator to determine the element id of the elemnt _point is in
+  virtual void updateCachedElementID();
+
+  /// number of PKAs to be generated per cm2-s
+  const unsigned int _flux;
+
+  /// time step
+  const Real _dt;
+
+  /// Mesh that comes from another generator
+  MooseMesh & _mesh;
+
+  /// boundary name
+  const BoundaryName _boundary;
+
+  /// surface area of boundary
+  const Real _boundary_surface_area;
+
+  /// the location from which to start PKAs
+  Point _point;
+
+  /// PKA nuclear charge
+  const unsigned int _Z;
+
+  /// PKA mass
+  const Real _m;
+
+  /// PKA Energy (in eV)
+  const Real _E;
+
+  /// cumulative probability list with element IDs
+  const std::vector<std::pair<Real, dof_id_type>> _prob_elem_pairs;
+
+  /// point locator to determine element pointers form locations
+  std::unique_ptr<PointLocatorBase> _pl;
+
+  /// the element id of the element containing _point
+  dof_id_type _elem_id;
+
+  Real boundarySurfaceArea(const BoundaryName & boundary);
+
+  std::vector<std::pair<Real, dof_id_type>> volumeWeightedElemDist(const BoundaryName & boundary);
+};

src/userobjects/PKASurfaceFluxGenerator.C

@@ -0,0 +1,184 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "PKASurfaceFluxGenerator.h"
+#include "MagpieUtils.h"
+
+registerMooseObject("MagpieApp", PKASurfaceFluxGenerator);
+
+InputParameters
+PKASurfaceFluxGenerator::validParams()
+{
+  InputParameters params = PKAGeneratorBase::validParams();
+  params.addClassDescription(
+      "This PKAGenerator starts particles from a random point in a boundary in a fixed direction.");
+  params.addRequiredParam<RealVectorValue>("direction", "The fixed direction the PKAs move along");
+  params.addRequiredParam<BoundaryName>("boundary", "The boundary to apply the flux to.");
+  params.addRequiredParam<unsigned int>(
+      "flux",
+      "The number of ions (starting PKAs) to strike the surface in ions/Area-second. Units of Area "
+      "must be consistent with either mesh unit length or given boundary_surface_area units.");
+  params.addRequiredParam<Real>("dt", "Time step size used in time stepper. Value in seconds.");
+  params.addParam<Real>(
+      "boundary_surface_area",
+      "Surface area of boundary. Used to determine number of initial "
+      "PKAs generated per timestep. Area must be consistent with flux dimensions.");
+  params.addRequiredParam<Real>("Z", "PKA nuclear charge");
+  params.addRequiredParam<Real>("m", "PKA mass in amu");
+  params.addRequiredParam<Real>("E", "PKA energy in eV");
+  return params;
+}
+
+PKASurfaceFluxGenerator::PKASurfaceFluxGenerator(const InputParameters & parameters)
+  : PKAGeneratorBase(parameters),
+    _direction(getParam<RealVectorValue>("direction")),
+    _flux(getParam<unsigned int>("flux")),
+    _dt(getParam<Real>("dt")),
+    _mesh(_fe_problem.mesh()),
+    _boundary(getParam<BoundaryName>("boundary")),
+    _boundary_surface_area(PKASurfaceFluxGenerator::boundarySurfaceArea(_boundary)),
+    _Z(getParam<Real>("Z")),
+    _m(getParam<Real>("m")),
+    _E(getParam<Real>("E")),
+    _prob_elem_pairs(PKASurfaceFluxGenerator::volumeWeightedElemDist(_boundary)),
+    _pl(_mesh.getPointLocator())
+{
+  _pl->enable_out_of_mesh_mode();
+  updateCachedElementID();
+}
+
+void
+PKASurfaceFluxGenerator::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_list,
+                                    const MyTRIMRasterizer::PKAParameters & pka_parameters,
+                                    const MyTRIMRasterizer::AveragedData &) const
+{
+  if (_current_elem->id() != _elem_id)
+    return;
+
+  unsigned int num_pka = _flux * _dt * _boundary_surface_area;
+  if (pka_parameters._recoil_rate_scaling != 1)
+    num_pka = std::floor(pka_parameters._recoil_rate_scaling * num_pka + getRandomReal());
+
+  for (unsigned i = 0; i < num_pka; ++i)
+  {
+    // each fission event generates a pair of recoils
+    MyTRIM_NS::IonBase pka;
+
+    // sample fission fragment masses
+    pka._Z = _Z;
+    pka._m = _m;
+    pka._E = _E;
+
+    // the tag is the element this PKA get registered as upon stopping
+    // -1 means the PKA will be ignored
+    pka._tag = ionTag(pka_parameters, pka._Z, pka._m);
+    ;
+
+    // set stopping criteria
+    pka.setEf();
+
+    // sample elements
+    Real rnd_num = getRandomReal();
+    auto it = std::lower_bound(_prob_elem_pairs.begin(),
+                               _prob_elem_pairs.end(),
+                               rnd_num,
+                               [](const std::pair<Real, dof_id_type> & entry, const Real & value)
+                               { return entry.first < value; });
+
+    // Get the selected element ID
+    dof_id_type selected_elem_id = it->second;
+
+    // Retrieve the element pointer from mesh
+    const Elem * rnd_elem = _mesh.getMesh().elem_ptr(selected_elem_id);
+
+    Point _point = MagpieUtils::randomElementPoint(
+        *rnd_elem, getRandomPoint()); // sample random point in element
+    pka._pos = _point;
+
+    // set random direction for ion 1 and opposite direction for ion 2
+    setDirection(pka);
+
+    // add PKA to list
+    ion_list.push_back(pka);
+  }
+}
+
+void
+PKASurfaceFluxGenerator::setDirection(MyTRIM_NS::IonBase & ion) const
+{
+  ion._dir = _direction;
+}
+
+Real
+PKASurfaceFluxGenerator::boundarySurfaceArea(const BoundaryName & boundary)
+{
+  if (isParamValid("boundary_surface_area"))
+    return getParam<Real>("boundary_surface_area");
+
+  BoundaryID boundary_id = _mesh.getBoundaryID(boundary);
+  Real volume_sum = 0;
+
+  const auto range = _mesh.getBoundaryElementRange();
+  for (const BndElement * bnd_elem : *range)
+  {
+    if (bnd_elem->_bnd_id != boundary_id)
+      continue; // skip other boundaries
+    const auto elem = bnd_elem->_elem;
+    const auto side = bnd_elem->_side;
+    volume_sum += elem->side_ptr(side)->volume();
+  }
+  mooseAssert(volume_sum > libMesh::TOLERANCE, "boundary_surface_area is not strictly positive!");
+  return volume_sum;
+}
+
+std::vector<std::pair<Real, dof_id_type>>
+PKASurfaceFluxGenerator::volumeWeightedElemDist(const BoundaryName & boundary)
+{
+  BoundaryID boundary_id = _mesh.getBoundaryID(boundary);
+  Real volume_sum = 0;
+  std::vector<std::pair<dof_id_type, Real>> elem_volumes; // (elem_id, side_volume)
+  std::vector<std::pair<Real, dof_id_type>> prob_elem_pairs;
+
+  const auto range = _mesh.getBoundaryElementRange();
+  for (const BndElement * bnd_elem : *range)
+  {
+    if (bnd_elem->_bnd_id != boundary_id)
+      continue; // skip other boundaries
+    const auto elem = bnd_elem->_elem;
+    dof_id_type elem_id = elem->id();
+    const auto side = bnd_elem->_side;
+    Real side_volume = elem->side_ptr(side)->volume();
+    volume_sum += side_volume;
+    elem_volumes.emplace_back(elem_id, side_volume);
+  }
+
+  // Normalize into cumulative distribution
+  Real cumulative_sum = 0.0;
+  for (const auto & [elem_id, side_vol] : elem_volumes)
+  {
+    cumulative_sum += side_vol / volume_sum;
+    prob_elem_pairs.emplace_back(cumulative_sum, elem_id);
+  }
+  mooseAssert(std::abs(cumulative_sum - 1.0) < libMesh::TOLERANCE,
+              "Volumes are not normalized to sum to 1.0!");
+  mooseAssert(std::abs(prob_elem_pairs.back().first - 1.0) < libMesh::TOLERANCE,
+              "Element probabilities are not normalized to sum to 1.0!");
+
+  return prob_elem_pairs;
+}
+
+void
+PKASurfaceFluxGenerator::updateCachedElementID()
+{
+  // get element containing the point
+  mooseAssert(_pl != nullptr, "initialize() must be called on the MooseMyTRIMSample object.");
+  const Elem * elem = (*_pl)(_point);
+  if (elem == nullptr)
+    mooseError("Point ", _point, " is not within the domain.");
+  _elem_id = elem->id();
+}

tests/userobjects/pka_surface_flux_generator/heat_deposition_2D.i

@@ -0,0 +1,123 @@
+#
+# ions depositing energy and heat conduction in 2D
+#
+
+[Mesh]
+  type = MyTRIMMesh
+  dim = 3
+  nx = 20
+  ny = 20
+  xmin = -100
+  xmax = 100
+  ymin = 0
+  ymax = 200
+[]
+
+[Variables]
+  [T]
+  []
+[]
+
+[Kernels]
+  [source]
+    type = MyTRIMElementHeatSource
+    variable = T
+    runner = runner
+  []
+  [conduction]
+    type = ADHeatConduction
+    variable = T
+    thermal_conductivity = 1 # dummy value
+  []
+  [time]
+    type = TimeDerivative
+    variable = T
+  []
+[]
+
+[AuxVariables]
+  [energy_dep]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+
+  [./cC]
+    order = CONSTANT
+    family = MONOMIAL
+    initial_condition = 1
+  [../]
+[]
+
+[AuxKernels]
+  [energy_dep]
+    variable = energy_dep
+    type = MyTRIMElementEnergyAux
+    runner = runner
+  []
+[]
+
+[UserObjects]
+  [./constant]
+    type = PKASurfaceFluxGenerator
+    E = 1e3
+    Z = 14
+    m = 28
+    direction = '1 0 0'
+    boundary = 'bottom'
+    flux = 1e3
+    dt = 1
+    boundary_surface_area = 1
+  [../]
+  [./rasterizer]
+    type = MyTRIMRasterizer
+    var = 'cC'
+    M = '12'
+    Z = '6'
+    Edisp = '16.3'
+    site_volume = 0.0404
+    pka_generator = constant
+    trim_module = ENERGY_DEPOSITION
+  [../]
+  [./runner]
+    type = MyTRIMElementRun
+    rasterizer = rasterizer
+  [../]
+[]
+
+[Postprocessors]
+  [E]
+    type = ElementIntegralVariablePostprocessor
+    variable = energy_dep
+    execute_on = timestep_end
+  []
+  [pka_total_E]
+    type = MyTRIMPKAInfo
+    rasterizer = rasterizer
+    value_type = TOTAL_ENERGY
+  []
+  [pka_total_Z]
+    type = MyTRIMPKAInfo
+    rasterizer = rasterizer
+    value_type = TOTAL_CHARGE
+  []
+  [pka_total_m]
+    type = MyTRIMPKAInfo
+    rasterizer = rasterizer
+    value_type = TOTAL_MASS
+  []
+  [pka_total_num]
+    type = MyTRIMPKAInfo
+    rasterizer = rasterizer
+    value_type = TOTAL_NUMBER
+  []
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 3
+  nl_abs_tol = 1e-10
+[]
+
+[Outputs]
+  exodus = true
+[]