Additional changes (idaholab#367)

include/userobjects/MDRunBase.h

@@ -113,6 +113,9 @@ class MDRunBase : public GeneralUserObject
   /// stores bounding boxes of all other processors
   std::vector<BoundingBox> _bbox;
 
+  /// maximum granular radius for parallel bounding boxes
+  Real _max_granular_radius;
+
   /// total number of particles
   unsigned int _n_particles = 0;
 
@@ -130,8 +133,4 @@ class MDRunBase : public GeneralUserObject
 
   /// A KDTree object to handle md_particles
   std::unique_ptr<KDTree> _kd_tree;
-
-  /// maximum granular radius
-  Real _max_granular_radius;
 };
-

src/userobjects/LAMMPSFileRunner.C

@@ -135,7 +135,10 @@ LAMMPSFileRunner::findBracketingLAMMPSFiles(Real md_time,
                  f.second,
                  " failed. LAMMPS filename must be base.<t>.xyz where <t> is the timestamp.");
     std::stringstream sstr(elements[1]);
-    sstr >> timestamp;
+
+    if (sstr >> timestamp || sstr.eof())
+      if (sstr.fail())
+        continue;
 
     // increase the counter & check if this file is a candidate for before/after
     ++lammps_files_found;
@@ -286,8 +289,12 @@ LAMMPSFileRunner::readLAMMPSFileHistory(std::pair<FileName, FileName> filenames,
   if (np != _n_particles)
     mooseError("Files ",
                filenames.first,
+               " : ",
+               _n_particles,
                " and ",
                filenames.second,
+               " : ",
+               np,
                " have different number of particles.");
 
   // skip another five lines

src/userobjects/MDRunBase.C

@@ -42,6 +42,11 @@ validParams<MDRunBase>()
   params.addParam<MultiMooseEnum>("md_particle_properties",
                                   MDRunBase::mdParticleProperties(),
                                   "Properties of MD particles to be obtained and stored.");
+  params.addRangeCheckedParam<Real>(
+      "max_granular_radius",
+      0,
+      "max_granular_radius>=0",
+      "Maximum radius of granular particles. This can be important for partitioning.");
   params.addClassDescription(
       "Base class for execution of coupled molecular dynamics MOOSE calculations.");
   return params;
@@ -56,6 +61,15 @@ MDRunBase::MDRunBase(const InputParameters & parameters)
     _nproc(_app.n_processors()),
     _bbox(_nproc)
 {
+  // if the calculation is granular, max_particle_radius
+  if (_granular)
+  {
+    if (!parameters.isParamSetByUser("max_granular_radius"))
+      paramError("max_granular_radius",
+                 "max_granular_radius must be set for granular calculations");
+    _max_granular_radius = getParam<Real>("max_granular_radius");
+  }
+
   // check _properties array and MooseEnum for consistency
   if (N_MD_PROPERTIES != mdParticleProperties().getNames().size())
     mooseError("Mismatch of MD particle property enum and property array. Report on github.");
@@ -68,22 +82,35 @@ MDRunBase::MDRunBase(const InputParameters & parameters)
 void
 MDRunBase::initialSetup()
 {
-  // TODO: need to inflate the bounding box for granular particles
-  if (_granular)
-    mooseDoOnce(mooseWarning("Granular particles can lead to wrong answers in parallel runs "
-                             "because processor bounding boxes do not account for particle size."));
-
   for (unsigned int j = 0; j < _nproc; ++j)
+  {
     _bbox[j] = MeshTools::create_processor_bounding_box(_mesh, j);
+
+    // inflate bounding box
+    for (unsigned int d = 0; d < _dim; ++d)
+    {
+      _bbox[j].first(d) -= _max_granular_radius;
+      _bbox[j].second(d) += _max_granular_radius;
+    }
+  }
 }
 
 void
 MDRunBase::timestepSetup()
 {
   // update/init subdomain bounding boxes
   for (unsigned int j = 0; j < _nproc; ++j)
+  {
     _bbox[j] = MeshTools::create_processor_bounding_box(_mesh, j);
 
+    // inflate bounding box
+    for (unsigned int d = 0; d < _dim; ++d)
+    {
+      _bbox[j].first(d) -= _max_granular_radius;
+      _bbox[j].second(d) += _max_granular_radius;
+    }
+  }
+
   // callback for updating md particle list
   updateParticleList();
 }
@@ -203,11 +230,27 @@ MDRunBase::updateElementGranularVolumes()
   // clear the granular volume map
   _elem_granular_volumes.clear();
 
-  // _max_granular_radius
+  /*
+     // Ideally we want to compute _max_granular_radius automatically
+     // but the value is needed for setting bounding boxes to partition
+     // the particles on the processors. For now, it's an input parameter
+     // until a path forward is determined.
+     // The commented lines compute the largest granular radius
   _max_granular_radius = 0;
   for (auto & p : _md_particles.properties)
     if (p[7] > _max_granular_radius)
       _max_granular_radius = p[7];
+  */
+  /// do a sanity check _max_granular_radius
+  Real mgr = 0;
+  for (auto & p : _md_particles.properties)
+    if (p[7] > mgr)
+      mgr = p[7];
+  if (mgr > _max_granular_radius)
+    mooseError("Granular particle with radius: ",
+               mgr,
+               " exceeds max_granular_radius: ",
+               _max_granular_radius);
 
   /// loop over all local elements
   ConstElemRange * active_local_elems = _mesh.getActiveLocalElementRange();

tests/userobjects/lammps_file/granular.i

@@ -53,6 +53,7 @@
     type = LAMMPSFileRunner
     lammps_file = 'granular.0.xyz'
     xyz_columns = '0 1 2'
+    max_granular_radius = 0.1
     md_particle_properties = 'radius charge'
     property_columns = '3 4'
   [../]

tests/userobjects/lammps_file/granular_sequence.i

@@ -44,6 +44,7 @@
     lammps_file = 'granular'
     time_sequence = true
     xyz_columns = '0 1 2'
+    max_granular_radius = 0.1
     md_particle_properties = 'radius charge'
     property_columns = '3 4'
   [../]

tests/userobjects/lammps_file/granular_tet.i

@@ -43,6 +43,7 @@
     type = LAMMPSFileRunner
     lammps_file = 'single.0.xyz'
     xyz_columns = '0 1 2'
+    max_granular_radius = 1    
     md_particle_properties = 'radius'
     property_columns = '3'
   [../]