Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
4 : /* */
5 : /* Copyright 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 :
9 : #include "LAMMPSFileRunner.h"
10 : #include "MooseUtils.h"
11 : #include "Function.h"
12 : #include <sstream>
13 : #include <cstring>
14 : #include <cctype>
15 :
16 : registerMooseObject("MagpieApp", LAMMPSFileRunner);
17 :
18 : InputParameters
19 48 : LAMMPSFileRunner::validParams()
20 : {
21 48 : InputParameters params = MDRunBase::validParams();
22 48 : params.addClassDescription("Reads a LAMMPS dump file or sequence of dump files and "
23 : "maps LAMMPS particles to MOOSE FEM mesh.");
24 96 : params.addParam<bool>("time_sequence",
25 96 : false,
26 : "If true, a sequence of dump files is read, "
27 : "else a single dump file is read.");
28 96 : params.addRequiredParam<FileName>("lammps_file",
29 : "Name of LAMMPS dump file or file base of LAMMPS file for "
30 : "a sequence.");
31 96 : params.addRequiredParam<std::vector<unsigned int>>(
32 : "xyz_columns", "Column ids of the x, y, and z coordinates of particles.");
33 96 : params.addParam<std::vector<unsigned int>>(
34 : "property_columns", {}, "Column ids of the properties.");
35 96 : params.addParam<FunctionName>("time_conversion",
36 : "A conversion from FEM simulation time to MD time stamps.");
37 48 : params.addClassDescription("Allows coupling FEM calculations to LAMMPS dump files.");
38 48 : return params;
39 0 : }
40 :
41 24 : LAMMPSFileRunner::LAMMPSFileRunner(const InputParameters & parameters)
42 : : MDRunBase(parameters),
43 24 : _time_sequence(getParam<bool>("time_sequence")),
44 24 : _lammps_file(getParam<FileName>("lammps_file")),
45 48 : _pos_columns(getParam<std::vector<unsigned int>>("xyz_columns")),
46 48 : _prop_columns(getParam<std::vector<unsigned int>>("property_columns")),
47 48 : _time_conversion(nullptr)
48 : {
49 : ;
50 24 : if (_pos_columns.size() != _dim)
51 0 : mooseError(
52 0 : "Dimensionality is ", _dim, " but xyz_columns has ", _pos_columns.size(), " entries.");
53 :
54 48 : if (isParamValid("time_conversion"))
55 8 : _time_conversion = &getFunction("time_conversion");
56 :
57 24 : if (_properties.size() != _prop_columns.size())
58 0 : mooseError("properties array must have same length as property_columns");
59 24 : }
60 :
61 : void
62 24 : LAMMPSFileRunner::initialSetup()
63 : {
64 : // call MDRunBase initialSetup
65 24 : MDRunBase::initialSetup();
66 :
67 24 : if (_time_sequence)
68 : return;
69 :
70 12 : readLAMMPSFile(_lammps_file);
71 12 : updateKDTree();
72 : }
73 :
74 : void
75 260 : LAMMPSFileRunner::updateParticleList()
76 : {
77 : // read LAMMPS file if a time sequence is read
78 260 : if (_time_sequence)
79 : {
80 : // map to MD time
81 248 : Real md_time = _t;
82 248 : if (_time_conversion)
83 240 : md_time = _time_conversion->value(_t, Point());
84 :
85 : // find files w/ timestamps "bracketing" md_time
86 : std::pair<std::string, std::string> files;
87 248 : std::pair<Real, Real> timestamps;
88 248 : findBracketingLAMMPSFiles(md_time, files, timestamps);
89 248 : readLAMMPSFileHistory(files, timestamps, md_time);
90 248 : updateKDTree();
91 : }
92 :
93 : // update the mapping to mesh if mesh has changed or ...
94 260 : mapMDParticles();
95 :
96 : // update the granular candidates as well if necessary
97 260 : if (_granular)
98 16 : updateElementGranularVolumes();
99 260 : }
100 :
101 : void
102 248 : LAMMPSFileRunner::findBracketingLAMMPSFiles(Real md_time,
103 : std::pair<std::string, std::string> & filenames,
104 : std::pair<Real, Real> & timestamps)
105 : {
106 : /// the return values
107 : std::string before;
108 : std::string after;
109 :
110 : /// extract path and file_base from provided lammps file entry
111 248 : std::pair<std::string, std::string> p = MooseUtils::splitFileName(_lammps_file);
112 :
113 : /// get all files in the directory
114 744 : std::list<std::string> files = MooseUtils::getFilesInDirs({p.first});
115 :
116 : /// loop over all files, retain the ones that begin with base and extract their timestamps
117 : std::vector<Real> times;
118 248 : times.reserve(files.size());
119 :
120 : unsigned int lammps_files_found = 0;
121 : Real closest_after_time = std::numeric_limits<Real>::max();
122 : Real closest_before_time = -std::numeric_limits<Real>::max();
123 1334 : for (auto & file : files)
124 : {
125 1086 : std::pair<std::string, std::string> f = MooseUtils::splitFileName(file);
126 1086 : if (f.second.compare(0, p.second.size(), p.second))
127 64 : continue;
128 :
129 : // split the file at dots and check the second item ==> timestamp
130 : Real timestamp;
131 : std::vector<std::string> elements;
132 2044 : MooseUtils::tokenize(f.second, elements, 1, ".");
133 1022 : if (elements.size() < 2)
134 0 : mooseError("Tokenizing filename ",
135 : f.second,
136 : " failed. LAMMPS filename must be base.<t>.xyz where <t> is the timestamp.");
137 1022 : std::stringstream sstr(elements[1]);
138 :
139 2044 : if (!isTimestamp(elements[1]))
140 : continue;
141 : sstr >> timestamp;
142 :
143 : // if (sstr >> timestamp || sstr.eof())
144 : // if (sstr.fail())
145 : // continue;
146 :
147 : // increase the counter & check if this file is a candidate for before/after
148 976 : ++lammps_files_found;
149 976 : if (timestamp > md_time && timestamp < closest_after_time)
150 : {
151 : closest_after_time = timestamp;
152 : after = file;
153 : }
154 668 : else if (timestamp <= md_time && timestamp > closest_before_time)
155 : {
156 : closest_before_time = timestamp;
157 : before = file;
158 : }
159 1022 : }
160 :
161 248 : if (lammps_files_found == 0)
162 0 : mooseError("No LAMMPS file with name base ", p.second, " in folder ", p.first);
163 :
164 : // guard against the case that only before or after were found => set the two equal
165 248 : if (before.size() == 0)
166 : before = after;
167 248 : else if (after.size() == 0)
168 : after = before;
169 :
170 496 : filenames = std::pair<std::string, std::string>(before, after);
171 : timestamps = std::pair<Real, Real>(closest_before_time, closest_after_time);
172 496 : }
173 :
174 : void
175 12 : LAMMPSFileRunner::readLAMMPSFile(FileName filename)
176 : {
177 : // check if this file is readable
178 12 : MooseUtils::checkFileReadable(filename);
179 :
180 : // read LAMMPS file
181 12 : std::ifstream file(filename);
182 :
183 : // skip first three lines
184 48 : for (unsigned int j = 0; j < 3; ++j)
185 36 : file.ignore(10000, '\n');
186 :
187 : // get number of particles
188 12 : file >> _n_particles;
189 :
190 : // skip another five lines
191 84 : for (unsigned int j = 0; j < 6; ++j)
192 72 : file.ignore(10000, '\n');
193 :
194 : // attach the right _md_particles.pos
195 12 : std::vector<Point> & position = _md_particles.pos;
196 :
197 : // zero out the position & id vector
198 : position.clear();
199 12 : _md_particles.id.clear();
200 12 : _md_particles.properties.clear();
201 :
202 : // size the md particle vector considering not all processors have all particles
203 12 : position.reserve(std::ceil(2.0 * _n_particles / _nproc));
204 12 : _md_particles.id.reserve(std::ceil(2.0 * _n_particles / _nproc));
205 12 : _md_particles.properties.reserve(std::ceil(2.0 * _n_particles / _nproc));
206 :
207 : // read particle entries
208 12 : _n_local_particles = 0;
209 56 : for (unsigned int j = 0; j < _n_particles; ++j)
210 : {
211 : // check if this particle is in the local bounding box
212 : std::string line;
213 44 : std::getline(file, line);
214 : std::vector<std::string> elements;
215 88 : MooseUtils::tokenize<std::string>(line, elements, 1, " ");
216 :
217 : #ifndef NDEBUG
218 : // check that enough columns exist in debug
219 : // find largest requested column
220 : unsigned int largest_column = 0;
221 : for (unsigned int i = 0; i < _dim; ++i)
222 : if (_pos_columns[i] > largest_column)
223 : largest_column = _pos_columns[i];
224 :
225 : for (unsigned int i = 0; i < _prop_columns.size(); ++i)
226 : if (_prop_columns[i] > largest_column)
227 : largest_column = _prop_columns[i];
228 :
229 : if (largest_column >= elements.size())
230 : mooseError("Error reading ", filename, " on line ", line);
231 : #endif
232 :
233 : // read location
234 : Point pos;
235 176 : for (unsigned int i = 0; i < _dim; ++i)
236 : {
237 132 : std::stringstream sstr(elements[_pos_columns[i]]);
238 : sstr >> pos(i);
239 132 : }
240 :
241 : // read properties
242 44 : std::vector<Real> props(_md_particles._prop_size);
243 88 : for (unsigned int i = 0; i < _prop_columns.size(); ++i)
244 : {
245 44 : unsigned int prop_id = _properties.get(i);
246 44 : std::stringstream sstr(elements[_prop_columns[i]]);
247 44 : sstr >> props[propIndex(prop_id)];
248 44 : }
249 :
250 : // check if this particle is in this processor BBox
251 44 : if (!_bbox.contains_point(pos))
252 : continue;
253 :
254 40 : ++_n_local_particles;
255 40 : position.push_back(pos);
256 40 : _md_particles.id.push_back(j);
257 40 : _md_particles.properties.push_back(props);
258 44 : }
259 12 : file.close();
260 :
261 : // size back
262 : position.shrink_to_fit();
263 : _md_particles.id.shrink_to_fit();
264 12 : }
265 :
266 : bool
267 1022 : LAMMPSFileRunner::isTimestamp(std::string ts_candidate) const
268 : {
269 2238 : for (auto & s : ts_candidate)
270 1262 : if (!isdigit(s))
271 : return false;
272 : return true;
273 : }
274 :
275 : void
276 248 : LAMMPSFileRunner::readLAMMPSFileHistory(std::pair<FileName, FileName> filenames,
277 : std::pair<Real, Real> timestamps,
278 : Real current_time)
279 : {
280 : // check if this file is readable
281 248 : MooseUtils::checkFileReadable(filenames.first);
282 248 : MooseUtils::checkFileReadable(filenames.second);
283 :
284 : // open LAMMPS file
285 248 : std::ifstream file_before(filenames.first);
286 248 : std::ifstream file_after(filenames.second);
287 :
288 : // compute weights for particle
289 248 : Real weight = (current_time - timestamps.first) / (timestamps.second - timestamps.first);
290 :
291 : // skip first three lines
292 992 : for (unsigned int j = 0; j < 3; ++j)
293 : {
294 744 : file_before.ignore(10000, '\n');
295 744 : file_after.ignore(10000, '\n');
296 : }
297 :
298 : // get number of particles
299 : unsigned int np;
300 248 : file_before >> _n_particles;
301 : file_after >> np;
302 248 : if (np != _n_particles)
303 0 : mooseError("Files ",
304 0 : filenames.first,
305 : " : ",
306 : _n_particles,
307 : " and ",
308 0 : filenames.second,
309 : " : ",
310 : np,
311 : " have different number of particles.");
312 :
313 : // skip another five lines
314 1736 : for (unsigned int j = 0; j < 6; ++j)
315 : {
316 1488 : file_before.ignore(10000, '\n');
317 1488 : file_after.ignore(10000, '\n');
318 : }
319 :
320 : // attach the right _md_particles.pos
321 248 : std::vector<Point> & position = _md_particles.pos;
322 :
323 : // zero out the position & id vector
324 : position.clear();
325 248 : _md_particles.id.clear();
326 248 : _md_particles.properties.clear();
327 :
328 : // size the md particle vector considering not all processors have all particles
329 248 : position.reserve(std::ceil(2.0 * _n_particles / _nproc));
330 248 : _md_particles.id.reserve(std::ceil(2.0 * _n_particles / _nproc));
331 248 : _md_particles.properties.reserve(std::ceil(2.0 * _n_particles / _nproc));
332 :
333 : // read particle entries
334 248 : _n_local_particles = 0;
335 648 : for (unsigned int j = 0; j < _n_particles; ++j)
336 : {
337 : // check if this particle is in the local bounding box
338 : // either in before or after state
339 :
340 : // get the MD particle from the before file
341 : std::vector<std::string> elements;
342 : std::string line_before;
343 400 : std::getline(file_before, line_before);
344 800 : MooseUtils::tokenize<std::string>(line_before, elements, 1, " ");
345 :
346 : #ifndef NDEBUG
347 : // check that enough columns exist in debug/devel
348 : // find largest requested column
349 : unsigned int largest_column = 0;
350 : for (unsigned int i = 0; i < _dim; ++i)
351 : if (_pos_columns[i] > largest_column)
352 : largest_column = _pos_columns[i];
353 :
354 : for (unsigned int i = 0; i < _prop_columns.size(); ++i)
355 : if (_prop_columns[i] > largest_column)
356 : largest_column = _prop_columns[i];
357 :
358 : if (largest_column >= elements.size())
359 : mooseError("Error reading ", filenames.first, " on line ", line_before);
360 : #endif
361 :
362 : Point pos_before;
363 1600 : for (unsigned int i = 0; i < _dim; ++i)
364 : {
365 1200 : std::stringstream sstr(elements[_pos_columns[i]]);
366 : sstr >> pos_before(i);
367 1200 : }
368 :
369 : // read properties from before file
370 400 : std::vector<Real> props_before(_md_particles._prop_size);
371 480 : for (unsigned int i = 0; i < _prop_columns.size(); ++i)
372 : {
373 80 : unsigned int prop_id = _properties.get(i);
374 80 : std::stringstream sstr(elements[_prop_columns[i]]);
375 80 : sstr >> props_before[propIndex(prop_id)];
376 80 : }
377 :
378 : // get the MD particle from the after file
379 : elements.clear();
380 : std::string line_after;
381 400 : std::getline(file_after, line_after);
382 :
383 800 : MooseUtils::tokenize<std::string>(line_after, elements, 1, " ");
384 :
385 : // we have determined largest_column in !NDEBUG so we can use it in mooseAssert
386 : mooseAssert(largest_column < elements.size(),
387 : "Error reading " << filenames.second << " on line " << line_after);
388 :
389 : Point pos_after;
390 1600 : for (unsigned int i = 0; i < _dim; ++i)
391 : {
392 1200 : std::stringstream sstr(elements[_pos_columns[i]]);
393 : sstr >> pos_after(i);
394 1200 : }
395 :
396 : // read properties from before file
397 400 : std::vector<Real> props_after(_md_particles._prop_size);
398 480 : for (unsigned int i = 0; i < _prop_columns.size(); ++i)
399 : {
400 80 : unsigned int prop_id = _properties.get(i);
401 80 : std::stringstream sstr(elements[_prop_columns[i]]);
402 80 : sstr >> props_after[propIndex(prop_id)];
403 80 : }
404 :
405 : // check if this particle is in this processor BBox
406 400 : if (!_bbox.contains_point(pos_after) && !_bbox.contains_point(pos_before))
407 : continue;
408 :
409 400 : ++_n_local_particles;
410 400 : position.push_back((1 - weight) * pos_before + weight * pos_after);
411 400 : _md_particles.id.push_back(j);
412 400 : _md_particles.properties.push_back(std::vector<Real>(_md_particles._prop_size));
413 480 : for (unsigned int i = 0; i < _md_particles._prop_size; ++i)
414 80 : _md_particles.properties.back()[i] = (1 - weight) * props_before[i] + weight * props_after[i];
415 400 : }
416 248 : file_before.close();
417 248 : file_after.close();
418 :
419 : // size back
420 : position.shrink_to_fit();
421 : _md_particles.id.shrink_to_fit();
422 248 : }
|
