1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing author: Norbert Podhorszki (Oak Ridge National Laboratory)
17 ------------------------------------------------------------------------- */
18
19 #include "reader_adios.h"
20 #include "comm.h"
21 #include "error.h"
22 #include "memory.h"
23 #include <cmath>
24 #include <cstdlib>
25 #include <cstring>
26
27 #include "adios2.h"
28 #include "math_const.h"
29
30 using namespace LAMMPS_NS;
31
32 using namespace MathConst;
33
34 // also in read_dump.cpp
35
36 enum { ID, TYPE, X, Y, Z, VX, VY, VZ, Q, IX, IY, IZ, FX, FY, FZ };
37 enum { UNSET, NOSCALE_NOWRAP, NOSCALE_WRAP, SCALE_NOWRAP, SCALE_WRAP };
38
39 #define SMALL 1.0e-6
40
41 // true if the difference between two floats is "small".
42 // cannot use fabsf() since it is not fully portable.
is_smalldiff(const float & val1,const float & val2)43 static bool is_smalldiff(const float &val1, const float &val2)
44 {
45 return (fabs(static_cast<double>(val1 - val2)) < SMALL);
46 }
47
48 namespace LAMMPS_NS
49 {
50 class ReadADIOSInternal
51 {
52
53 public:
ReadADIOSInternal()54 ReadADIOSInternal() {};
55 ~ReadADIOSInternal() = default;
56
57 // name of adios group, referrable in adios2_config.xml
58 const std::string ioName = "read_dump";
59 adios2::ADIOS *ad = nullptr; // adios object
60 adios2::IO io; // adios group of variables and attributes in this dump
61 adios2::Engine fh; // adios file/stream handle object
62 // ADIOS input variables we need to change every step
63 adios2::Variable<uint64_t> varNtimestep;
64 adios2::Variable<uint64_t> varNatoms;
65 adios2::Variable<double> varAtoms;
66 // list of column names for the atom table
67 // (individual list of 'columns' string)
68 std::vector<std::string> columnNames;
69 float timeout = 0.0;
70 };
71 } // namespace LAMMPS_NS
72
73 /* ---------------------------------------------------------------------- */
74
ReaderADIOS(LAMMPS * lmp)75 ReaderADIOS::ReaderADIOS(LAMMPS *lmp) : Reader(lmp)
76 {
77 fieldindex = nullptr;
78 nAtoms = 0;
79 nAtomsTotal = 0;
80 atomOffset = 0;
81 nstep = 0;
82 nid = 0;
83 me = comm->me;
84
85 internal = new ReadADIOSInternal();
86 try {
87 internal->ad =
88 new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
89 } catch (std::ios_base::failure &e) {
90 char str[256];
91 snprintf(str, sizeof(str), "ADIOS initialization failed with error: %s",
92 e.what());
93 error->one(FLERR, str);
94 }
95
96 /* Define the group holding all variables and attributes */
97 internal->io = internal->ad->DeclareIO(internal->ioName);
98 }
99
100 /* ---------------------------------------------------------------------- */
101
~ReaderADIOS()102 ReaderADIOS::~ReaderADIOS()
103 {
104 if (me == 0) {
105 memory->destroy(fieldindex);
106 }
107 internal->columnNames.clear();
108 if (internal->fh) {
109 internal->fh.Close();
110 }
111 delete internal->ad;
112 delete internal;
113 }
114
115 /* ----------------------------------------------------------------------
116 pass on settings to find and load the proper plugin
117 Called by all processors.
118 ------------------------------------------------------------------------- */
settings(int narg,char ** arg)119 void ReaderADIOS::settings(int narg, char **arg)
120 {
121 int idx = 0;
122 while (idx < narg) {
123 if (!strcmp(arg[idx], "timeout")) {
124 if (idx + 1 < narg) {
125 internal->timeout = std::stof(arg[idx + 1]);
126 internal->io.SetParameter("OpenTimeoutSecs", arg[idx + 1]);
127 ++idx;
128 } else {
129 char str[128];
130 snprintf(str, sizeof(str),
131 "Missing value for 'timeout' option for ADIOS "
132 "read_dump command");
133 error->one(FLERR, str);
134 }
135 }
136 ++idx;
137 }
138 }
139
140 /* ----------------------------------------------------------------------
141 try to open given file
142 Every process must call this Collective operation
143 ------------------------------------------------------------------------- */
144
open_file(const char * file)145 void ReaderADIOS::open_file(const char *file)
146 {
147 int rv;
148 char str[1024];
149
150 // close open file, if needed.
151 if (internal->fh)
152 internal->fh.Close();
153
154 try {
155 internal->fh = internal->io.Open(file, adios2::Mode::Read, world);
156 } catch (std::ios_base::failure &e) {
157 char str[256];
158 snprintf(str, sizeof(str), "%s", e.what());
159 error->one(FLERR, str);
160 }
161 if (!internal->fh) {
162 snprintf(str, sizeof(str), "Cannot open file %s using ADIOS", file);
163 error->one(FLERR, str);
164 }
165 }
166
167 /* ----------------------------------------------------------------------
168 close current file
169 Every process must call this Collective operation
170 ------------------------------------------------------------------------- */
171
close_file()172 void ReaderADIOS::close_file()
173 {
174 // close open file, if needed.
175 if (internal->fh) {
176 internal->fh.Close();
177 }
178 }
179
180 /* ----------------------------------------------------------------------
181 read and return time stamp from dump file
182 if first read reaches end-of-file, return 1 so caller can open next file
183 Called by all processors.
184 ------------------------------------------------------------------------- */
185
read_time(bigint & ntimestep)186 int ReaderADIOS::read_time(bigint &ntimestep)
187 {
188 char str[1024];
189
190 adios2::StepStatus status =
191 internal->fh.BeginStep(adios2::StepMode::Read, internal->timeout);
192
193 switch (status) {
194 case adios2::StepStatus::EndOfStream:
195 case adios2::StepStatus::NotReady:
196 case adios2::StepStatus::OtherError:
197 return 1;
198 default:
199 break;
200 }
201
202 internal->varNtimestep =
203 internal->io.InquireVariable<uint64_t>("ntimestep");
204
205 if (!internal->varNtimestep) {
206 snprintf(str, sizeof(str),
207 "Did not find 'ntimestep' variable in ADIOS file %s",
208 internal->fh.Name().c_str());
209 error->one(FLERR, str);
210 }
211
212 ntimestep = static_cast<bigint>(internal->varNtimestep.Max());
213 // std::cerr << " **** ReaderADIOS::read_time found step " << ntimestep
214 // << " **** " << std::endl;
215 return 0;
216 }
217
218 /* ----------------------------------------------------------------------
219 skip snapshot for given timestamp
220 Called by all processors.
221 ------------------------------------------------------------------------- */
222
skip()223 void ReaderADIOS::skip() { internal->fh.EndStep(); }
224
225 /* ----------------------------------------------------------------------
226 read remaining header info:
227 return natoms
228 box bounds, triclinic (inferred), fieldflag (1 if any fields not found),
229 xyz flags = from input scaleflag & wrapflag
230 if fieldflag set:
231 match Nfield fields to per-atom column labels
232 allocate and set fieldindex = which column each field maps to
233 fieldtype = X,VX,IZ etc
234 fieldlabel = user-specified label or nullptr if use fieldtype default
235 xyz flag = scaledflag if has fieldlabel name, else set by x,xs,xu,xsu
236 only called by proc 0
237 ------------------------------------------------------------------------- */
238
read_header(double box[3][3],int & boxinfo,int & triclinic,int fieldinfo,int nfield,int * fieldtype,char ** fieldlabel,int scaleflag,int wrapflag,int & fieldflag,int & xflag,int & yflag,int & zflag)239 bigint ReaderADIOS::read_header(double box[3][3], int &boxinfo, int &triclinic,
240 int fieldinfo, int nfield, int *fieldtype,
241 char **fieldlabel, int scaleflag, int wrapflag,
242 int &fieldflag, int &xflag, int &yflag,
243 int &zflag)
244 {
245 char str[1024];
246 nid = 0;
247
248 // signal that we have no box info at all so far.
249
250 internal->varNatoms = internal->io.InquireVariable<uint64_t>("natoms");
251 if (!internal->varNatoms) {
252 snprintf(str, sizeof(str),
253 "Did not find 'natoms' variable in ADIOS file %s",
254 internal->fh.Name().c_str());
255 error->one(FLERR, str);
256 }
257
258 /* nAtoms */
259 nAtomsTotal = internal->varNatoms.Max();
260 uint64_t rem = nAtomsTotal % comm->nprocs;
261 nAtoms = nAtomsTotal / comm->nprocs;
262 atomOffset = comm->me * nAtoms;
263 if (comm->me < rem) {
264 ++nAtoms;
265 atomOffset += comm->me;
266 } else {
267 atomOffset += rem;
268 }
269
270 /* triclinic */
271 adios2::Attribute<int32_t> attTriclinic =
272 internal->io.InquireAttribute<int32_t>("triclinic");
273 if (!attTriclinic) {
274 snprintf(str, sizeof(str),
275 "Did not find 'triclinic' attribute in ADIOS file %s",
276 internal->fh.Name().c_str());
277 error->one(FLERR, str);
278 }
279
280 triclinic = attTriclinic.Data()[0];
281
282 /* read Box */
283 adios2::Variable<double> varBoxxlo =
284 internal->io.InquireVariable<double>("boxxlo");
285 adios2::Variable<double> varBoxxhi =
286 internal->io.InquireVariable<double>("boxxhi");
287 adios2::Variable<double> varBoxylo =
288 internal->io.InquireVariable<double>("boxylo");
289 adios2::Variable<double> varBoxyhi =
290 internal->io.InquireVariable<double>("boxyhi");
291 adios2::Variable<double> varBoxzlo =
292 internal->io.InquireVariable<double>("boxzlo");
293 adios2::Variable<double> varBoxzhi =
294 internal->io.InquireVariable<double>("boxzhi");
295
296 box[0][0] = varBoxxlo.Max();
297 box[0][1] = varBoxxhi.Max();
298 box[0][2] = 0.0;
299 box[1][0] = varBoxylo.Max();
300 box[1][1] = varBoxyhi.Max();
301 box[1][2] = 0.0;
302 box[2][0] = varBoxzlo.Max();
303 box[2][1] = varBoxzhi.Max();
304 box[2][2] = 0.0;
305
306 if (triclinic) {
307 adios2::Variable<double> varBoxxy =
308 internal->io.InquireVariable<double>("boxxy");
309 adios2::Variable<double> varBoxxz =
310 internal->io.InquireVariable<double>("boxxz");
311 adios2::Variable<double> varBoxyz =
312 internal->io.InquireVariable<double>("boxyz");
313
314 box[0][2] = varBoxxy.Max();
315 box[1][2] = varBoxxz.Max();
316 box[2][2] = varBoxyz.Max();
317 }
318
319 boxinfo = 1;
320
321 // if no field info requested, just return
322
323 if (!fieldinfo)
324 return nAtoms;
325
326 memory->create(fieldindex, nfield, "read_dump:fieldindex");
327
328 /* Columns */
329 adios2::Attribute<std::string> attColumns =
330 internal->io.InquireAttribute<std::string>("columns");
331
332 std::vector<std::string> labelVector = attColumns.Data();
333 int nwords = labelVector.size();
334 std::map<std::string, int> labels;
335 for (int i = 0; i < nwords; ++i) {
336 labels.emplace(labelVector[i], i);
337 }
338
339 int s_index, u_index, su_index;
340 xflag = UNSET;
341 yflag = UNSET;
342 zflag = UNSET;
343
344 // copy fieldtype list for supported fields
345
346 for (int i = 0; i < nfield; i++) {
347 if (fieldlabel[i]) {
348 fieldindex[i] = find_label(fieldlabel[i], labels);
349 if (fieldtype[i] == X)
350 xflag = 2 * scaleflag + wrapflag + 1;
351 else if (fieldtype[i] == Y)
352 yflag = 2 * scaleflag + wrapflag + 1;
353 else if (fieldtype[i] == Z)
354 zflag = 2 * scaleflag + wrapflag + 1;
355 }
356
357 else if (fieldtype[i] == ID)
358 fieldindex[i] = find_label("id", labels);
359 else if (fieldtype[i] == TYPE)
360 fieldindex[i] = find_label("type", labels);
361
362 else if (fieldtype[i] == X) {
363 fieldindex[i] = find_label("x", labels);
364 xflag = NOSCALE_WRAP;
365 if (fieldindex[i] < 0) {
366 fieldindex[i] = nwords;
367 s_index = find_label("xs", labels);
368 u_index = find_label("xu", labels);
369 su_index = find_label("xsu", labels);
370 if (s_index >= 0 && s_index < fieldindex[i]) {
371 fieldindex[i] = s_index;
372 xflag = SCALE_WRAP;
373 }
374 if (u_index >= 0 && u_index < fieldindex[i]) {
375 fieldindex[i] = u_index;
376 xflag = NOSCALE_NOWRAP;
377 }
378 if (su_index >= 0 && su_index < fieldindex[i]) {
379 fieldindex[i] = su_index;
380 xflag = SCALE_NOWRAP;
381 }
382 }
383 if (fieldindex[i] == nwords)
384 fieldindex[i] = -1;
385
386 } else if (fieldtype[i] == Y) {
387 fieldindex[i] = find_label("y", labels);
388 yflag = NOSCALE_WRAP;
389 if (fieldindex[i] < 0) {
390 fieldindex[i] = nwords;
391 s_index = find_label("ys", labels);
392 u_index = find_label("yu", labels);
393 su_index = find_label("ysu", labels);
394 if (s_index >= 0 && s_index < fieldindex[i]) {
395 fieldindex[i] = s_index;
396 yflag = SCALE_WRAP;
397 }
398 if (u_index >= 0 && u_index < fieldindex[i]) {
399 fieldindex[i] = u_index;
400 yflag = NOSCALE_NOWRAP;
401 }
402 if (su_index >= 0 && su_index < fieldindex[i]) {
403 fieldindex[i] = su_index;
404 yflag = SCALE_NOWRAP;
405 }
406 }
407 if (fieldindex[i] == nwords)
408 fieldindex[i] = -1;
409
410 } else if (fieldtype[i] == Z) {
411 fieldindex[i] = find_label("z", labels);
412 zflag = NOSCALE_WRAP;
413 if (fieldindex[i] < 0) {
414 fieldindex[i] = nwords;
415 s_index = find_label("zs", labels);
416 u_index = find_label("zu", labels);
417 su_index = find_label("zsu", labels);
418 if (s_index >= 0 && s_index < fieldindex[i]) {
419 fieldindex[i] = s_index;
420 zflag = SCALE_WRAP;
421 }
422 if (u_index >= 0 && u_index < fieldindex[i]) {
423 fieldindex[i] = u_index;
424 zflag = NOSCALE_NOWRAP;
425 }
426 if (su_index >= 0 && su_index < fieldindex[i]) {
427 fieldindex[i] = su_index;
428 zflag = SCALE_NOWRAP;
429 }
430 }
431 if (fieldindex[i] == nwords)
432 fieldindex[i] = -1;
433
434 } else if (fieldtype[i] == VX)
435 fieldindex[i] = find_label("vx", labels);
436 else if (fieldtype[i] == VY)
437 fieldindex[i] = find_label("vy", labels);
438 else if (fieldtype[i] == VZ)
439 fieldindex[i] = find_label("vz", labels);
440
441 else if (fieldtype[i] == FX)
442 fieldindex[i] = find_label("fx", labels);
443 else if (fieldtype[i] == FY)
444 fieldindex[i] = find_label("fy", labels);
445 else if (fieldtype[i] == FZ)
446 fieldindex[i] = find_label("fz", labels);
447
448 else if (fieldtype[i] == Q)
449 fieldindex[i] = find_label("q", labels);
450
451 else if (fieldtype[i] == IX)
452 fieldindex[i] = find_label("ix", labels);
453 else if (fieldtype[i] == IY)
454 fieldindex[i] = find_label("iy", labels);
455 else if (fieldtype[i] == IZ)
456 fieldindex[i] = find_label("iz", labels);
457 }
458
459 // set fieldflag = -1 if any unfound fields
460
461 fieldflag = 0;
462 for (int i = 0; i < nfield; i++)
463 if (fieldindex[i] < 0)
464 fieldflag = -1;
465
466 return nAtoms;
467 }
468
469 /* ----------------------------------------------------------------------
470 read N atom lines from dump file
471 stores appropriate values in fields array
472 return 0 if success, 1 if error
473 only called by proc 0
474 ------------------------------------------------------------------------- */
475
read_atoms(int n,int nfield,double ** fields)476 void ReaderADIOS::read_atoms(int n, int nfield, double **fields)
477 {
478 char str[1024];
479
480 /* Read Atoms */
481 /* This is the firsts (and last) read of array data, so we
482 * call EndStep() here instead of PerformGets()
483 */
484
485 adios2::Variable<double> varAtoms =
486 internal->io.InquireVariable<double>("atoms");
487
488 if (n != nAtoms) {
489 snprintf(
490 str, sizeof(str),
491 "ReaderADIOS::read_atoms() expects 'n=%d' equal to the number of "
492 "atoms (=%" PRIu64 ") for process %d in ADIOS file %s.",
493 n, nAtoms, comm->me, internal->fh.Name().c_str());
494 error->one(FLERR, str);
495 }
496
497 size_t ncols = varAtoms.Count()[1];
498 varAtoms.SetSelection({{atomOffset, 0}, {nAtoms, ncols}});
499
500 std::vector<double> table;
501
502 internal->fh.Get<double>(varAtoms, table);
503 // EndStep or PerformGets required to make the read happen
504 internal->fh.EndStep();
505
506 size_t idx;
507 for (int i = 0; i < nAtoms; i++) {
508 idx = i * ncols;
509 for (int m = 0; m < nfield; m++) {
510 fields[i][m] = table[idx + fieldindex[m]];
511 }
512 }
513 }
514
515 /* ----------------------------------------------------------------------
516 match label to any of N labels
517 return index of match or -1 if no match
518 ------------------------------------------------------------------------- */
519
find_label(const std::string & label,const std::map<std::string,int> & labels)520 int ReaderADIOS::find_label(const std::string &label,
521 const std::map<std::string, int> &labels)
522 {
523 std::map<std::string, int>::const_iterator it = labels.find(label);
524 if (it != labels.end()) {
525 return it->second;
526 }
527 return -1;
528 }
529