First design of MD capability in Magpie (idaholab#367)

include/userobjects/LAMMPSFileRunner.h

@@ -0,0 +1,54 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#ifndef LAMMPSFILERUNNER_H
+#define LAMMPSFILERUNNER_H
+
+#include "MDRunBase.h"
+
+class LAMMPSFileRunner;
+
+template <>
+InputParameters validParams<LAMMPSFileRunner>();
+
+/**
+ * Reads lammps dump files to emulate MD simulation
+ */
+class LAMMPSFileRunner : public MDRunBase
+{
+public:
+  LAMMPSFileRunner(const InputParameters & parameters);
+
+  virtual void initialize() override {}
+  virtual void execute() override {}
+  virtual void finalize() override {}
+  virtual void initialSetup() override;
+
+  virtual void updateParticleList() override;
+
+protected:
+  /// reads a LAMMPS file into pos or pos_old
+  void readLAMMPSFile(FileName filename, bool read_old_pos = false);
+
+  /// whether a sequence of input files or a single input file is read
+  bool _time_sequence;
+
+  /// name of LAMMPS file or file base if _time_sequence == true
+  FileName _lammps_file;
+
+  /// column of , y, z coordinate in LAMMPS files
+  std::array<unsigned int, 3> _pos_columns;
+
+  /// get velocity from file or compute it from positions
+  bool _velocity_from_file;
+
+  /// column of , y, z coordinate in LAMMPS files
+  std::array<unsigned int, 3> _vel_columns;
+};
+
+#endif // LAMMPSFILERUNNER_H

include/userobjects/MDRunBase.h

@@ -0,0 +1,91 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#ifndef MDRUNBASE_H
+#define MDRUNBASE_H
+
+#include "GeneralUserObject.h"
+#include "KDTree.h"
+#include "libmesh/bounding_box.h"
+
+class MDRunBase;
+class MooseMesh;
+
+template <>
+InputParameters validParams<MDRunBase>();
+
+/**
+ * Base class for molecular dynamics runs in Magpie
+ */
+class MDRunBase : public GeneralUserObject
+{
+public:
+  MDRunBase(const InputParameters & parameters);
+
+  void initialSetup() override;
+  void timestepSetup() override;
+
+  class MDParticles
+  {
+  public:
+    ///@{ particle's position and history
+    std::vector<Point> pos;
+    std::vector<Point> pos_old;
+    ///@}
+
+    ///@{ particle's velocity and history
+    std::vector<Point> vel;
+    std::vector<Point> vel_old;
+    ///@}
+
+    /// the id of particle in the MD calculation
+    std::vector<unsigned int> id;
+
+    /// the mesh element the particles are in
+    std::vector<dof_id_type> elem_id;
+  };
+
+protected:
+  /// call back function to update the particle list
+  virtual void updateParticleList() = 0;
+
+  /// updates the KDTree object
+  void updateKDTree();
+
+  /// map MDParticles to elements
+  void mapMDParticles();
+
+  /// The Mesh we're using
+  MooseMesh & _mesh;
+
+  /// dimension of the mesh
+  const unsigned int _dim;
+
+  /// dimension of the mesh
+  const unsigned int _nproc;
+
+  /// stores bounding boxes of all other processors
+  std::vector<BoundingBox> _bbox;
+
+  /// total number of particles
+  unsigned int _n_particles;
+
+  /// number of local particles
+  unsigned int _n_local_particles;
+
+  /// stores the location of
+  MDParticles _md_particles;
+
+  /// a map from elem pointer to particles in this element
+  std::map<dof_id_type, std::vector<unsigned int>> _elem_particles;
+
+  /// A KDTree object to handle md_particles
+  std::unique_ptr<KDTree> _kd_tree;
+};
+
+#endif // MDRUNBASE_H

src/userobjects/LAMMPSFileRunner.C

@@ -0,0 +1,133 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "LAMMPSFileRunner.h"
+#include "MooseUtils.h"
+#include <sstream>
+
+registerMooseObject("MagpieApp", LAMMPSFileRunner);
+
+template <>
+InputParameters
+validParams<LAMMPSFileRunner>()
+{
+  InputParameters params = validParams<MDRunBase>();
+  params.addClassDescription("Reads a LAMMPS dump file or sequence of dump files and "
+                             "maps LAMMPS particles to MOOSE FEM mesh.");
+  params.addParam<bool>("time_sequence",
+                        false,
+                        "If true, a sequence of dump files is read, "
+                        "else a single dump file is read.");
+  params.addRequiredParam<FileName>("lammps_file",
+                                    "Name of LAMMPS dump file or file base of LAMMPS file for "
+                                    "a sequence.");
+  params.addRequiredParam<std::vector<unsigned int>>(
+      "xyz_columns", "Column ids of the x, y, and z coordinates of particles.");
+
+  return params;
+}
+
+LAMMPSFileRunner::LAMMPSFileRunner(const InputParameters & parameters)
+  : MDRunBase(parameters),
+    _time_sequence(getParam<bool>("time_sequence")),
+    _lammps_file(getParam<FileName>("lammps_file"))
+{
+  std::vector<unsigned int> pos = getParam<std::vector<unsigned int>>("xyz_columns");
+  if (pos.size() != _dim)
+    mooseError("Dimensionality is ", _dim, " but xyz_columns has ", pos.size(), " entries.");
+  for (unsigned int j = 0; j < _dim; ++j)
+    _pos_columns[j] = pos[j];
+}
+
+void
+LAMMPSFileRunner::initialSetup()
+{
+  // call MDRunBase initialSetup
+  MDRunBase::initialSetup();
+
+  // read the file only if we don't read a sequence of files
+  if (_time_sequence)
+    return;
+
+  readLAMMPSFile(_lammps_file);
+  updateKDTree();
+}
+
+void
+LAMMPSFileRunner::updateParticleList()
+{
+  // read LAMMPS file if a time sequence is read
+  if (!_time_sequence)
+  {
+  }
+
+  // update the mapping to mesh if mesh has changed or ...
+  mapMDParticles();
+
+  unsigned int np = 0;
+  for (auto & p : _elem_particles)
+    np += p.second.size();
+}
+
+void
+LAMMPSFileRunner::readLAMMPSFile(FileName filename, bool read_old_pos)
+{
+  // read LAMMPS file
+  std::ifstream file(filename);
+
+  // skip first three lines
+  for (unsigned int j = 0; j < 3; ++j)
+    file.ignore(10000, '\n');
+
+  // get number of particles
+  file >> _n_particles;
+
+  // skip another five lines
+  for (unsigned int j = 0; j < 6; ++j)
+    file.ignore(10000, '\n');
+
+  // attach the right _md_particles.pos
+  std::vector<Point> & position = _md_particles.pos;
+  if (read_old_pos)
+    position = _md_particles.pos_old;
+
+  // size the md particle vector considering not all processors have all particles
+  position.reserve(std::ceil(2.0 * _n_particles / _nproc));
+  _md_particles.id.reserve(std::ceil(2.0 * _n_particles / _nproc));
+
+  // read particle entries
+  _n_local_particles = 0;
+  for (unsigned int j = 0; j < _n_particles; ++j)
+  {
+    // check if this particle is in the local bounding box
+    std::string line;
+    std::getline(file, line);
+    std::vector<std::string> elements;
+    MooseUtils::tokenize<std::string>(line, elements, 1, " ");
+
+    Point pos;
+    for (unsigned int i = 0; i < _dim; ++i)
+    {
+      std::stringstream sstr(elements[_pos_columns[i]]);
+      sstr >> pos(i);
+    }
+
+    // check if this particle is in this processor BBox
+    if (!_bbox[processor_id()].contains_point(pos))
+      continue;
+
+    ++_n_local_particles;
+    position.push_back(pos);
+    _md_particles.id.push_back(j);
+  }
+  file.close();
+
+  // size back
+  position.shrink_to_fit();
+  _md_particles.id.shrink_to_fit();
+}