Implement read_dump with format 'adios'. It modifies read_dump.cpp to support... (580909fd) · Commits · 郑智淋 / lammps

examples/USER/adios/in.adios_balance

+3 −2

Original line number	Diff line number	Diff line
		@@ -39,8 +39,9 @@ fix 10 all balance 50 0.9 rcb

		compute 1 all property/atom proc
		variable p atom c_1%10
		dump 2 all custom 50 balance.dump id v_p x y z
		dump 3 all custom/adios 50 balance_custom.bp id v_p x y z
		dump 5 all atom 50 balance_atom.dump
		dump 2 all custom 50 balance_custom.dump id type x y z vx vy vz ix iy iz fx fy fz
		dump 3 all custom/adios 50 balance_custom.bp id type x y z vx vy vz ix iy iz fx fy fz
		dump 4 all atom/adios 50 balance_atom.bp

		#dump 3 all image 50 image.*.jpg v_p type &

examples/USER/adios/in.adios_balance2

0 → 100644

+60 −0

Original line number	Diff line number	Diff line
		# 2d circle of particles inside a box with LJ walls

		variable b index 0

		variable x index 50
		variable y index 20
		variable d index 20
		variable v index 5
		variable w index 2

		units lj
		dimension 2
		atom_style atomic
		boundary f f p

		lattice hex 0.85
		region box block 0 $x 0 $y -0.5 0.5
		create_box 1 box
		region circle sphere $(v_d/2+1) $(v_d/2/sqrt(3.0)+1) 0.0 $(v_d/2)
		create_atoms 1 region circle
		mass 1 1.0

		velocity all create 0.5 87287 loop geom
		velocity all set $v $w 0 sum yes

		pair_style lj/cut 2.5
		pair_coeff 1 1 10.0 1.0 2.5

		neighbor 0.3 bin
		neigh_modify delay 0 every 1 check yes

		fix 1 all nve

		fix 2 all wall/lj93 xlo 0.0 1 1 2.5 xhi $x 1 1 2.5
		fix 3 all wall/lj93 ylo 0.0 1 1 2.5 yhi $y 1 1 2.5

		comm_style tiled
		fix 10 all balance 50 0.9 rcb

		#read_dump balance_custom.dump 200 x y
		read_dump balance_atom.bp 200 x y format adios
		compute 1 all property/atom proc
		variable p atom c_1%10
		dump 2 all custom 50 balance2.dump id v_p x y z
		dump 3 all custom/adios 50 balance_custom2.bp id v_p x y z
		dump 4 all atom/adios 50 balance_atom2.bp

		#dump 3 all image 50 image.*.jpg v_p type &
		# adiam 1.0 view 0 0 zoom 1.8 subbox yes 0.02
		#variable colors string &
		# "red green blue yellow white &
		# purple pink orange lime gray"
		#dump_modify 3 pad 5 amap 0 10 sa 1 10 ${colors}

		thermo_style custom step temp epair press f_10[3] f_10
		thermo 100

		run 200
		write_dump all atom/adios balance_atom_final2.bp
		write_dump all atom balance_atom_final2.dump

src/USER-ADIOS/dump_atom_adios.cpp

+17 −8

Original line number	Diff line number	Diff line
		@@ -16,7 +16,6 @@
		------------------------------------------------------------------------- */

		#include "dump_atom_adios.h"
		#include <cstring>
		#include "atom.h"
		#include "domain.h"
		#include "error.h"
		@@ -24,6 +23,7 @@
		#include "memory.h"
		#include "universe.h"
		#include "update.h"
		#include <cstring>

		#include "adios2.h"

		@@ -50,7 +50,7 @@ public:
		// one ADIOS output variable we need to change every step
		adios2::Variable<double> varAtoms;
		};
		}
		} // namespace LAMMPS_NS

		/* ---------------------------------------------------------------------- */

		@@ -248,14 +248,21 @@ void DumpAtomADIOS::init_style()

		// setup column string

		if (scale_flag == 0 && image_flag == 0)
		std::vector<std::string> columnNames;

		if (scale_flag == 0 && image_flag == 0) {
		columns = (char *)"id type x y z";
		else if (scale_flag == 0 && image_flag == 1)
		columnNames = {"id", "type", "x", "y", "z"};
		} else if (scale_flag == 0 && image_flag == 1) {
		columns = (char *)"id type x y z ix iy iz";
		else if (scale_flag == 1 && image_flag == 0)
		columnNames = {"id", "type", "x", "y", "z", "ix", "iy", "iz"};
		} else if (scale_flag == 1 && image_flag == 0) {
		columns = (char *)"id type xs ys zs";
		else if (scale_flag == 1 && image_flag == 1)
		columnNames = {"id", "type", "xs", "ys", "zs"};
		} else if (scale_flag == 1 && image_flag == 1) {
		columns = (char *)"id type xs ys zs ix iy iz";
		columnNames = {"id", "type", "xs", "ys", "zs", "ix", "iy", "iz"};
		}

		// setup function ptrs

		@@ -316,7 +323,10 @@ void DumpAtomADIOS::init_style()
		int boundaryptr = reinterpret_cast<int >(domain->boundary);
		internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);

		internal->io.DefineAttribute<std::string>("columns", columns);
		size_t nColumns = static_cast<size_t>(size_one);
		internal->io.DefineAttribute<std::string>("columns", columnNames.data(),
		nColumns);
		internal->io.DefineAttribute<std::string>("columnstr", columns);
		internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
		internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "atom");
		internal->io.DefineAttribute<std::string>("LAMMPS/version",
		@@ -332,7 +342,6 @@ void DumpAtomADIOS::init_style()
		// atom table size is not known at the moment
		// it will be correctly defined at the moment of write
		size_t UnknownSizeYet = 1;
		size_t nColumns = static_cast<size_t>(size_one);
		internal->varAtoms = internal->io.DefineVariable<double>(
		"atoms", {UnknownSizeYet, nColumns}, {UnknownSizeYet, 0},
		{UnknownSizeYet, nColumns});

src/USER-ADIOS/dump_custom_adios.cpp

+4 −4

Original line number	Diff line number	Diff line
		@@ -16,8 +16,6 @@
		------------------------------------------------------------------------- */

		#include "dump_custom_adios.h"
		#include <cmath>
		#include <cstring>
		#include "atom.h"
		#include "compute.h"
		#include "domain.h"
		@@ -32,6 +30,8 @@
		#include "universe.h"
		#include "update.h"
		#include "variable.h"
		#include <cmath>
		#include <cstring>

		#include "adios2.h"

		@@ -124,7 +124,7 @@ public:
		// (individual list of 'columns' string)
		std::vector<std::string> columnNames;
		};
		}
		} // namespace LAMMPS_NS

		/* ---------------------------------------------------------------------- */

		@@ -415,7 +415,7 @@ void DumpCustomADIOS::init_style()
		"columns", internal->columnNames.data(), nColumns);
		internal->io.DefineAttribute<std::string>("columnstr", columns);
		internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
		internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "atom");
		internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "custom");
		internal->io.DefineAttribute<std::string>("LAMMPS/version",
		universe->version);
		internal->io.DefineAttribute<std::string>("LAMMPS/num_ver",

src/USER-ADIOS/reader_adios.cpp

0 → 100644

+491 −0

Original line number	Diff line number	Diff line
		/* ----------------------------------------------------------------------
		LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
		http://lammps.sandia.gov, Sandia National Laboratories
		Steve Plimpton, sjplimp@sandia.gov

		Copyright (2003) Sandia Corporation. Under the terms of Contract
		DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
		certain rights in this software. This software is distributed under
		the GNU General Public License.

		See the README file in the top-level LAMMPS directory.
		------------------------------------------------------------------------- */

		/* ----------------------------------------------------------------------
		Contributing author: Axel Kohlmeyer (Temple U)
		------------------------------------------------------------------------- */

		#include "reader_adios.h"
		#include "comm.h"
		#include "error.h"
		#include "memory.h"
		#include <cmath>
		#include <cstdlib>

		#include "adios2.h"
		#include "math_const.h"

		using namespace LAMMPS_NS;

		using namespace MathConst;

		// also in read_dump.cpp

		enum { ID, TYPE, X, Y, Z, VX, VY, VZ, Q, IX, IY, IZ, FX, FY, FZ };
		enum { UNSET, NOSCALE_NOWRAP, NOSCALE_WRAP, SCALE_NOWRAP, SCALE_WRAP };

		#define SMALL 1.0e-6

		// true if the difference between two floats is "small".
		// cannot use fabsf() since it is not fully portable.
		static bool is_smalldiff(const float &val1, const float &val2)
		{
		return (fabs(static_cast<double>(val1 - val2)) < SMALL);
		}

		namespace LAMMPS_NS
		{
		class ReadADIOSInternal
		{

		public:
		ReadADIOSInternal(){};
		~ReadADIOSInternal() = default;

		// name of adios group, referrable in adios2_config.xml
		const std::string ioName = "read_dump";
		adios2::ADIOS *ad = nullptr; // adios object
		adios2::IO io; // adios group of variables and attributes in this dump
		adios2::Engine fh; // adios file/stream handle object
		// ADIOS input variables we need to change every step
		adios2::Variable<uint64_t> varNtimestep;
		adios2::Variable<uint64_t> varNatoms;
		adios2::Variable<double> varAtoms;
		// list of column names for the atom table
		// (individual list of 'columns' string)
		std::vector<std::string> columnNames;
		};
		} // namespace LAMMPS_NS

		/* ---------------------------------------------------------------------- */

		ReaderADIOS::ReaderADIOS(LAMMPS *lmp) : Reader(lmp)
		{
		fieldindex = NULL;
		nAtoms = 0;
		nAtomsTotal = 0;
		atomOffset = 0;
		nstep = 0;
		nid = 0;
		me = comm->me;

		internal = new ReadADIOSInternal();
		internal->ad =
		new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);

		/* Define the group holding all variables and attributes */
		internal->io = internal->ad->DeclareIO(internal->ioName);
		}

		/* ---------------------------------------------------------------------- */

		ReaderADIOS::~ReaderADIOS()
		{
		if (me == 0) {
		memory->destroy(fieldindex);
		}
		internal->columnNames.clear();
		if (internal->fh) {
		internal->fh.Close();
		}
		delete internal->ad;
		delete internal;
		}

		/* ----------------------------------------------------------------------
		pass on settings to find and load the proper plugin
		Called by all processors.
		------------------------------------------------------------------------- */
		void ReaderADIOS::settings(int narg, char **arg) {}

		/* ----------------------------------------------------------------------
		try to open given file
		Every process must call this Collective operation
		------------------------------------------------------------------------- */

		void ReaderADIOS::open_file(const char *file)
		{
		int rv;
		char str[1024];

		// close open file, if needed.
		if (internal->fh)
		internal->fh.Close();

		internal->fh = internal->io.Open(file, adios2::Mode::Read, world);
		if (!internal->fh) {
		snprintf(str, strlen(str), "Cannot open file %s using ADIOS", file);
		error->one(FLERR, str);
		}
		}

		/* ----------------------------------------------------------------------
		close current file
		Every process must call this Collective operation
		------------------------------------------------------------------------- */

		void ReaderADIOS::close_file()
		{
		// close open file, if needed.
		if (internal->fh) {
		internal->fh.Close();
		}
		}

		/* ----------------------------------------------------------------------
		read and return time stamp from dump file
		if first read reaches end-of-file, return 1 so caller can open next file
		Called by all processors.
		------------------------------------------------------------------------- */

		int ReaderADIOS::read_time(bigint &ntimestep)
		{
		char str[1024];

		adios2::StepStatus status =
		internal->fh.BeginStep(adios2::StepMode::Read, 10.0f);

		switch (status) {
		case adios2::StepStatus::EndOfStream:
		case adios2::StepStatus::NotReady:
		case adios2::StepStatus::OtherError:
		return 1;
		default:
		break;
		}

		internal->varNtimestep =
		internal->io.InquireVariable<uint64_t>("ntimestep");

		if (!internal->varNtimestep) {
		snprintf(str, strlen(str),
		"Did not find 'ntimestep' variable in ADIOS file %s",
		internal->fh.Name().c_str());
		error->one(FLERR, str);
		}

		ntimestep = static_cast<bigint>(internal->varNtimestep.Max());
		return 0;
		}

		/* ----------------------------------------------------------------------
		skip snapshot for given timestamp
		Called by all processors.
		------------------------------------------------------------------------- */

		void ReaderADIOS::skip() { internal->fh.EndStep(); }

		/* ----------------------------------------------------------------------
		read remaining header info:
		return natoms
		box bounds, triclinic (inferred), fieldflag (1 if any fields not found),
		xyz flags = from input scaleflag & wrapflag
		if fieldflag set:
		match Nfield fields to per-atom column labels
		allocate and set fieldindex = which column each field maps to
		fieldtype = X,VX,IZ etc
		fieldlabel = user-specified label or NULL if use fieldtype default
		xyz flag = scaledflag if has fieldlabel name, else set by x,xs,xu,xsu
		only called by proc 0
		------------------------------------------------------------------------- */

		bigint ReaderADIOS::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)
		{
		char str[1024];
		nid = 0;

		// signal that we have no box info at all so far.

		internal->varNatoms = internal->io.InquireVariable<uint64_t>("natoms");
		if (!internal->varNatoms) {
		snprintf(str, strlen(str),
		"Did not find 'natoms' variable in ADIOS file %s",
		internal->fh.Name().c_str());
		error->one(FLERR, str);
		}

		/* nAtoms */
		nAtomsTotal = internal->varNatoms.Max();
		uint64_t rem = nAtomsTotal % comm->nprocs;
		nAtoms = nAtomsTotal / comm->nprocs;
		atomOffset = comm->me * nAtoms;
		if (comm->me < rem) {
		++nAtoms;
		atomOffset += comm->me;
		} else {
		atomOffset += rem;
		}

		/* triclinic */
		adios2::Attribute<int32_t> attTriclinic =
		internal->io.InquireAttribute<int32_t>("triclinic");
		if (!attTriclinic) {
		snprintf(str, strlen(str),
		"Did not find 'triclinic' attribute in ADIOS file %s",
		internal->fh.Name().c_str());
		error->one(FLERR, str);
		}

		triclinic = attTriclinic.Data()[0];

		/* read Box */
		adios2::Variable<double> varBoxxlo =
		internal->io.InquireVariable<double>("boxxlo");
		adios2::Variable<double> varBoxxhi =
		internal->io.InquireVariable<double>("boxxhi");
		adios2::Variable<double> varBoxylo =
		internal->io.InquireVariable<double>("boxylo");
		adios2::Variable<double> varBoxyhi =
		internal->io.InquireVariable<double>("boxyhi");
		adios2::Variable<double> varBoxzlo =
		internal->io.InquireVariable<double>("boxzlo");
		adios2::Variable<double> varBoxzhi =
		internal->io.InquireVariable<double>("boxzhi");

		box[0][0] = varBoxxlo.Max();
		box[0][1] = varBoxxhi.Max();
		box[0][2] = 0.0;
		box[1][0] = varBoxylo.Max();
		box[1][1] = varBoxyhi.Max();
		box[1][2] = 0.0;
		box[2][0] = varBoxzlo.Max();
		box[2][1] = varBoxzhi.Max();
		box[2][2] = 0.0;

		if (triclinic) {
		adios2::Variable<double> varBoxxy =
		internal->io.InquireVariable<double>("boxxy");
		adios2::Variable<double> varBoxxz =
		internal->io.InquireVariable<double>("boxxz");
		adios2::Variable<double> varBoxyz =
		internal->io.InquireVariable<double>("boxyz");

		box[0][2] = varBoxxy.Max();
		box[1][2] = varBoxxz.Max();
		box[2][2] = varBoxyz.Max();
		}

		boxinfo = 1;

		// if no field info requested, just return

		if (!fieldinfo)
		return nAtoms;

		memory->create(fieldindex, nfield, "read_dump:fieldindex");

		/* Columns */
		adios2::Attribute<std::string> attColumns =
		internal->io.InquireAttribute<std::string>("columns");

		std::vector<std::string> labelVector = attColumns.Data();
		int nwords = labelVector.size();
		std::map<std::string, int> labels;
		for (int i = 0; i < nwords; ++i) {
		labels.emplace(labelVector[i], i);
		}

		int s_index, u_index, su_index;
		xflag = UNSET;
		yflag = UNSET;
		zflag = UNSET;

		// copy fieldtype list for supported fields

		for (int i = 0; i < nfield; i++) {
		if (fieldlabel[i]) {
		fieldindex[i] = find_label(fieldlabel[i], labels);
		if (fieldtype[i] == X)
		xflag = 2 * scaleflag + wrapflag + 1;
		else if (fieldtype[i] == Y)
		yflag = 2 * scaleflag + wrapflag + 1;
		else if (fieldtype[i] == Z)
		zflag = 2 * scaleflag + wrapflag + 1;
		}

		else if (fieldtype[i] == ID)
		fieldindex[i] = find_label("id", labels);
		else if (fieldtype[i] == TYPE)
		fieldindex[i] = find_label("type", labels);

		else if (fieldtype[i] == X) {
		fieldindex[i] = find_label("x", labels);
		xflag = NOSCALE_WRAP;
		if (fieldindex[i] < 0) {
		fieldindex[i] = nwords;
		s_index = find_label("xs", labels);
		u_index = find_label("xu", labels);
		su_index = find_label("xsu", labels);
		if (s_index >= 0 && s_index < fieldindex[i]) {
		fieldindex[i] = s_index;
		xflag = SCALE_WRAP;
		}
		if (u_index >= 0 && u_index < fieldindex[i]) {
		fieldindex[i] = u_index;
		xflag = NOSCALE_NOWRAP;
		}
		if (su_index >= 0 && su_index < fieldindex[i]) {
		fieldindex[i] = su_index;
		xflag = SCALE_NOWRAP;
		}
		}
		if (fieldindex[i] == nwords)
		fieldindex[i] = -1;

		} else if (fieldtype[i] == Y) {
		fieldindex[i] = find_label("y", labels);
		yflag = NOSCALE_WRAP;
		if (fieldindex[i] < 0) {
		fieldindex[i] = nwords;
		s_index = find_label("ys", labels);
		u_index = find_label("yu", labels);
		su_index = find_label("ysu", labels);
		if (s_index >= 0 && s_index < fieldindex[i]) {
		fieldindex[i] = s_index;
		yflag = SCALE_WRAP;
		}
		if (u_index >= 0 && u_index < fieldindex[i]) {
		fieldindex[i] = u_index;
		yflag = NOSCALE_NOWRAP;
		}
		if (su_index >= 0 && su_index < fieldindex[i]) {
		fieldindex[i] = su_index;
		yflag = SCALE_NOWRAP;
		}
		}
		if (fieldindex[i] == nwords)
		fieldindex[i] = -1;

		} else if (fieldtype[i] == Z) {
		fieldindex[i] = find_label("z", labels);
		zflag = NOSCALE_WRAP;
		if (fieldindex[i] < 0) {
		fieldindex[i] = nwords;
		s_index = find_label("zs", labels);
		u_index = find_label("zu", labels);
		su_index = find_label("zsu", labels);
		if (s_index >= 0 && s_index < fieldindex[i]) {
		fieldindex[i] = s_index;
		zflag = SCALE_WRAP;
		}
		if (u_index >= 0 && u_index < fieldindex[i]) {
		fieldindex[i] = u_index;
		zflag = NOSCALE_NOWRAP;
		}
		if (su_index >= 0 && su_index < fieldindex[i]) {
		fieldindex[i] = su_index;
		zflag = SCALE_NOWRAP;
		}
		}
		if (fieldindex[i] == nwords)
		fieldindex[i] = -1;

		} else if (fieldtype[i] == VX)
		fieldindex[i] = find_label("vx", labels);
		else if (fieldtype[i] == VY)
		fieldindex[i] = find_label("vy", labels);
		else if (fieldtype[i] == VZ)
		fieldindex[i] = find_label("vz", labels);

		else if (fieldtype[i] == FX)
		fieldindex[i] = find_label("fx", labels);
		else if (fieldtype[i] == FY)
		fieldindex[i] = find_label("fy", labels);
		else if (fieldtype[i] == FZ)
		fieldindex[i] = find_label("fz", labels);

		else if (fieldtype[i] == Q)
		fieldindex[i] = find_label("q", labels);

		else if (fieldtype[i] == IX)
		fieldindex[i] = find_label("ix", labels);
		else if (fieldtype[i] == IY)
		fieldindex[i] = find_label("iy", labels);
		else if (fieldtype[i] == IZ)
		fieldindex[i] = find_label("iz", labels);
		}

		// set fieldflag = -1 if any unfound fields

		fieldflag = 0;
		for (int i = 0; i < nfield; i++)
		if (fieldindex[i] < 0)
		fieldflag = -1;

		return nAtoms;
		}

		/* ----------------------------------------------------------------------
		read N atom lines from dump file
		stores appropriate values in fields array
		return 0 if success, 1 if error
		only called by proc 0
		------------------------------------------------------------------------- */

		void ReaderADIOS::read_atoms(int n, int nfield, double **fields)
		{
		char str[1024];

		/* Read Atoms */
		/* This is the firsts (and last) read of array data, so we
		* call EndStep() here instead of PerformGets()
		*/

		adios2::Variable<double> varAtoms =
		internal->io.InquireVariable<double>("atoms");

		if (n != nAtoms) {
		snprintf(
		str, strlen(str),
		"ReaderADIOS::read_atoms() expects 'n=%d' equal to the number of "
		"atoms (=%" PRIu64 ") for process %d in ADIOS file %s.",
		n, nAtoms, comm->me, internal->fh.Name().c_str());
		error->one(FLERR, str);
		}

		size_t ncols = varAtoms.Count()[1];
		varAtoms.SetSelection({{atomOffset, 0}, {nAtoms, ncols}});

		std::vector<double> table;

		internal->fh.Get<double>(varAtoms, table);
		// EndStep or PerformGets required to make the read happen
		internal->fh.EndStep();

		size_t idx;
		for (int i = 0; i < nAtoms; i++) {
		idx = i * ncols;
		for (int m = 0; m < nfield; m++) {
		fields[i][m] = table[idx + fieldindex[m]];
		}
		}
		}

		/* ----------------------------------------------------------------------
		match label to any of N labels
		return index of match or -1 if no match
		------------------------------------------------------------------------- */

		int ReaderADIOS::find_label(const std::string &label,
		const std::map<std::string, int> &labels)
		{
		std::map<std::string, int>::const_iterator it = labels.find(label);
		if (it != labels.end()) {
		return it->second;
		}
		return -1;
		}

Admin message