diff --git a/include/userobjects/MDRunBase.h b/include/userobjects/MDRunBase.h index 4ab5e9d8..7cd619a4 100644 --- a/include/userobjects/MDRunBase.h +++ b/include/userobjects/MDRunBase.h @@ -113,6 +113,9 @@ class MDRunBase : public GeneralUserObject /// stores bounding boxes of all other processors std::vector _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 _kd_tree; - - /// maximum granular radius - Real _max_granular_radius; }; - diff --git a/src/userobjects/LAMMPSFileRunner.C b/src/userobjects/LAMMPSFileRunner.C index d709eb4c..cd7cff0b 100644 --- a/src/userobjects/LAMMPSFileRunner.C +++ b/src/userobjects/LAMMPSFileRunner.C @@ -135,7 +135,10 @@ LAMMPSFileRunner::findBracketingLAMMPSFiles(Real md_time, f.second, " failed. LAMMPS filename must be base..xyz where 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 filenames, if (np != _n_particles) mooseError("Files ", filenames.first, + " : ", + _n_particles, " and ", filenames.second, + " : ", + np, " have different number of particles."); // skip another five lines diff --git a/src/userobjects/MDRunBase.C b/src/userobjects/MDRunBase.C index 901c7028..0e074010 100644 --- a/src/userobjects/MDRunBase.C +++ b/src/userobjects/MDRunBase.C @@ -42,6 +42,11 @@ validParams() params.addParam("md_particle_properties", MDRunBase::mdParticleProperties(), "Properties of MD particles to be obtained and stored."); + params.addRangeCheckedParam( + "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("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,13 +82,17 @@ 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 @@ -82,8 +100,17 @@ 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(); diff --git a/tests/userobjects/lammps_file/granular.i b/tests/userobjects/lammps_file/granular.i index beebb2f7..0347a4d3 100644 --- a/tests/userobjects/lammps_file/granular.i +++ b/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' [../] diff --git a/tests/userobjects/lammps_file/granular_sequence.i b/tests/userobjects/lammps_file/granular_sequence.i index 2ff0eeda..eb8afb32 100644 --- a/tests/userobjects/lammps_file/granular_sequence.i +++ b/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' [../] diff --git a/tests/userobjects/lammps_file/granular_tet.i b/tests/userobjects/lammps_file/granular_tet.i index e471de6f..eddb2193 100644 --- a/tests/userobjects/lammps_file/granular_tet.i +++ b/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' [../]