Implement energy deposition tallying in Dirac runs (#396)

include/other/ThreadedRecoilDiracSourceLoop.h

@@ -39,10 +39,7 @@ class ThreadedRecoilDiracSourceLoop : public ThreadedRecoilLoopBase
   void addDefectToResult(const Point & p, unsigned int var, Real weight, DefectType type);
 
   /// add deposited energy to the result list
-  void addEnergyToResult(const Point & /*p*/, Real /*edep*/)
-  {
-    mooseError("Energy deposition is not implemented for dirac-based TRIM runs.");
-  }
+  void addEnergyToResult(const Point & p, Real edep);
 
   /// data such as interstitials and vacancies produced will be stored here
   MyTRIMResultList _result_list;
@@ -53,22 +50,23 @@ struct ThreadedRecoilDiracSourceLoop::MyTRIMResult
   MyTRIMResult(const Point & location,
                unsigned int var,
                ThreadedRecoilDiracSourceLoop::DefectType type,
-               dof_id_type elem_id)
-    : _location(location), _var(var), _type(type), _elem_id(elem_id)
+               dof_id_type elem_id,
+               Real weight = 1)
+    : _location(location), _var(var), _type(type), _elem_id(elem_id), _weight(weight)
   {
   }
 
-  MyTRIMResult() : _location(), _var(0), _type(NONE), _elem_id(libMesh::invalid_uint) {}
+  MyTRIMResult() : _location(), _var(0), _type(NONE), _elem_id(libMesh::invalid_uint), _weight(1) {}
 
   Point _location;
   unsigned int _var;
   DefectType _type;
   dof_id_type _elem_id;
+  Real _weight;
 };
 
 template <>
 void dataStore(std::ostream &, ThreadedRecoilDiracSourceLoop::MyTRIMResult &, void *);
 
 template <>
 void dataLoad(std::istream &, ThreadedRecoilDiracSourceLoop::MyTRIMResult &, void *);
-

include/other/ThreadedRecoilLoopBase.h

@@ -61,15 +61,16 @@ class ThreadedRecoilLoopBase
                         Real & distance) const;
 #endif // GSL_ENABLED
 
-  /// defect type enum (vacancies, interstitials, and incoming and outgoing replacements)
+  /// defect type enum (vacancies, interstitials, incoming and outgoing replacements, and energy deposition)
   enum DefectType
   {
     NONE = -1,
     VACANCY = 0,
     INTERSTITIAL,
     REPLACEMENT_IN,
     REPLACEMENT_OUT,
-    N_DEFECTS
+    N_DEFECTS,
+    ENERGY_DEPOSITION = 5
   };
 
 protected:

include/vectorpostprocessors/MyTRIMDiracEnergyResult.h

@@ -0,0 +1,42 @@
+/**********************************************************************/
+/*                     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 "GeneralVectorPostprocessor.h"
+
+class MyTRIMDiracEnergyResult;
+class MyTRIMDiracRun;
+
+template <>
+InputParameters validParams<MyTRIMDiracEnergyResult>();
+
+/**
+ * Outputs the list of MyTRIM energy deposition events
+ * compiled with the the MyTRIMDiracRunner
+ */
+class MyTRIMDiracEnergyResult : public GeneralVectorPostprocessor
+{
+public:
+  MyTRIMDiracEnergyResult(const InputParameters & parameters);
+
+  virtual void initialize();
+  virtual void execute();
+  virtual void finalize();
+
+protected:
+  const MyTRIMDiracRun & _mytrim;
+
+  ///@{ coordinates of the defects
+  VectorPostprocessorValue & _x;
+  VectorPostprocessorValue & _y;
+  VectorPostprocessorValue & _z;
+  VectorPostprocessorValue & _elem_id;
+  VectorPostprocessorValue & _energy_deposition;
+  ///@}
+};

src/other/ThreadedRecoilDiracSourceLoop.C

@@ -43,6 +43,14 @@ ThreadedRecoilDiracSourceLoop::addDefectToResult(const Point & p,
     _result_list.push_back(MyTRIMResult(p, var, type, elem->id()));
 }
 
+void
+ThreadedRecoilDiracSourceLoop::addEnergyToResult(const Point & p, Real edep)
+{
+  const Elem * elem = (*_pl)(p);
+  if (elem != nullptr)
+    _result_list.push_back(MyTRIMResult(p, 0, ENERGY_DEPOSITION, elem->id(), edep));
+}
+
 template <>
 void
 dataStore(std::ostream & stream, ThreadedRecoilDiracSourceLoop::MyTRIMResult & dsl, void * context)
@@ -51,6 +59,7 @@ dataStore(std::ostream & stream, ThreadedRecoilDiracSourceLoop::MyTRIMResult & d
   dataStore(stream, dsl._var, context);
   dataStore(stream, dsl._type, context);
   dataStore(stream, dsl._elem_id, context);
+  dataStore(stream, dsl._weight, context);
 }
 
 template <>
@@ -61,4 +70,5 @@ dataLoad(std::istream & stream, ThreadedRecoilDiracSourceLoop::MyTRIMResult & ds
   dataLoad(stream, dsl._var, context);
   dataLoad(stream, dsl._type, context);
   dataLoad(stream, dsl._elem_id, context);
+  dataLoad(stream, dsl._weight, context);
 }

src/vectorpostprocessors/MyTRIMDiracEnergyResult.C

@@ -0,0 +1,63 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "MyTRIMDiracEnergyResult.h"
+#include "MyTRIMDiracRun.h"
+
+registerMooseObject("MagpieApp", MyTRIMDiracEnergyResult);
+
+template <>
+InputParameters
+validParams<MyTRIMDiracEnergyResult>()
+{
+  InputParameters params = validParams<GeneralVectorPostprocessor>();
+  params.addRequiredParam<UserObjectName>(
+      "runner", "Name of the MyTRIMDiracRun userobject to pull data from.");
+  return params;
+}
+
+MyTRIMDiracEnergyResult::MyTRIMDiracEnergyResult(const InputParameters & parameters)
+  : GeneralVectorPostprocessor(parameters),
+    _mytrim(getUserObject<MyTRIMDiracRun>("runner")),
+    _x(declareVector("x")),
+    _y(declareVector("y")),
+    _z(declareVector("z")),
+    _elem_id(declareVector("elem_id")),
+    _energy_deposition(declareVector("energy_deposition"))
+{
+}
+
+void
+MyTRIMDiracEnergyResult::initialize()
+{
+  _x.clear();
+  _y.clear();
+  _z.clear();
+  _elem_id.clear();
+  _energy_deposition.clear();
+}
+
+void
+MyTRIMDiracEnergyResult::execute()
+{
+  for (auto && entry : _mytrim.result())
+    if (static_cast<unsigned int>(entry._type) == 5)
+    {
+      _x.push_back(entry._location(0));
+      _y.push_back(entry._location(1));
+      _z.push_back(entry._location(2));
+      _elem_id.push_back(entry._elem_id);
+      _energy_deposition.push_back(entry._weight);
+    }
+}
+
+void
+MyTRIMDiracEnergyResult::finalize()
+{
+  // the MyTRIMDiracRunner already does the necessary communication
+}

tests/vectorpostprocessors/dirac_energy_result.i

@@ -0,0 +1,76 @@
+box_half_length = 50
+
+[Mesh]
+  type = MyTRIMMesh
+  dim = 3
+  nx = 10
+  ny = 10
+  nz = 10
+  xmin = -${box_half_length}
+  xmax =  ${box_half_length}
+  ymin = -${box_half_length}
+  ymax =  ${box_half_length}
+  zmin = -${box_half_length}
+  zmax =  ${box_half_length}
+[]
+
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [./T]
+  [../]
+[]
+
+[AuxVariables]
+  # unit is atoms per nm^3
+  [./c]
+    initial_condition = 84.91232
+  [../]
+[]
+
+[UserObjects]
+  [./pka_gun]
+    type = PKAGun
+    m = 63.55
+    Z = 29
+    E = 50e3
+    direction = '1 0 0'
+    point = '0 0 0'
+    num_pkas = 1
+  [../]
+  [./rasterizer]
+    type = MyTRIMRasterizer
+    var = c
+    M = '63.55'
+    Z = '29'
+    pka_generator = pka_gun
+    length_unit = NANOMETER
+    var_physical_meaning = NUMBER_DENSITY
+    trim_module = ENERGY_DEPOSITION
+  [../]
+  [./runner]
+    type = MyTRIMDiracRun
+    rasterizer = rasterizer
+  [../]
+[]
+
+[VectorPostprocessors]
+  [./dirac_energy]
+    type = MyTRIMDiracEnergyResult
+    runner = runner
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 1
+  nl_abs_tol = 1e-10
+[]
+
+[Outputs]
+  exodus = true
+  csv = true
+  hide = c
+[]