Merge pull request #1296 from Adrian-Diaz/multi-comm-tiled (787d7d28) · Commits · 郑智淋 / lammps

doc/src/pair_colloid.rst

+1 −1

Original line number	Diff line number	Diff line
		@@ -77,7 +77,7 @@ The colloid-solvent interaction energy is given by
		\left[ 1 - \frac{\left(5 ~ a^6+45~a^4~r^2+63~a^2~r^4+15~r^6\right) \sigma^6}
		{15 \left(a-r\right)^6 \left( a+r \right)^6} \right], \quad r < r_c

		where :math:A_{cs}` is the Hamaker constant, a is the radius of the colloidal
		where :math:`A_{cs}` is the Hamaker constant, a is the radius of the colloidal
		particle, and :math:`r_c` is the cutoff. This formula is derived from the
		colloid-colloid interaction, letting one of the particle sizes go to
		zero.

src/atom_vec.cpp

+4 −2

Original line number	Diff line number	Diff line
		@@ -1096,8 +1096,9 @@ void AtomVec::unpack_border(int n, int first, double *buf)

		m = 0;
		last = first + n;
		while (last > nmax) grow(0);

		for (i = first; i < last; i++) {
		if (i == nmax) grow(0);
		x[i][0] = buf[m++];
		x[i][1] = buf[m++];
		x[i][2] = buf[m++];
		@@ -1165,8 +1166,9 @@ void AtomVec::unpack_border_vel(int n, int first, double *buf)

		m = 0;
		last = first + n;
		while (last > nmax) grow(0);

		for (i = first; i < last; i++) {
		if (i == nmax) grow(0);
		x[i][0] = buf[m++];
		x[i][1] = buf[m++];
		x[i][2] = buf[m++];

src/comm_tiled.cpp

+287 −114

Original line number	Diff line number	Diff line
		@@ -11,6 +11,11 @@
		See the README file in the top-level LAMMPS directory.
		------------------------------------------------------------------------- */

		/* ----------------------------------------------------------------------
		Contributing author (multi and triclinic support):
		Adrian Diaz (University of Florida)
		------------------------------------------------------------------------- */

		#include "comm_tiled.h"
		#include <mpi.h>
		#include <cmath>
		@@ -42,6 +47,11 @@ CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)
		style = 1;
		layout = Comm::LAYOUT_UNIFORM;
		pbc_flag = NULL;
		buf_send = NULL;
		buf_recv = NULL;
		overlap = NULL;
		rcbinfo = NULL;
		cutghostmulti = NULL;
		init_buffers();
		}

		@@ -69,6 +79,7 @@ CommTiled::~CommTiled()
		memory->destroy(overlap);
		deallocate_swap(maxswap);
		memory->sfree(rcbinfo);
		memory->destroy(cutghostmulti);
		}

		/* ----------------------------------------------------------------------
		@@ -84,11 +95,12 @@ void CommTiled::init_buffers()

		maxoverlap = 0;
		overlap = NULL;
		rcbinfo = NULL;
		cutghostmulti = NULL;
		sendbox_multi = NULL;

		maxswap = 6;
		allocate_swap(maxswap);

		rcbinfo = NULL;
		}

		/* ---------------------------------------------------------------------- */
		@@ -97,16 +109,18 @@ void CommTiled::init()
		{
		Comm::init();

		// cannot set nswap in init_buffers() b/c
		// dimension command can be after comm_style command

		nswap = 2*domain->dimension;

		memory->destroy(cutghostmulti);
		if (mode == Comm::MULTI)
		memory->create(cutghostmulti,atom->ntypes+1,3,"comm:cutghostmulti");

		int bufextra_old = bufextra;
		init_exchange();
		if (bufextra > bufextra_old) grow_send(maxsend+bufextra,2);

		// temporary restrictions

		if (mode == Comm::MULTI)
		error->all(FLERR,"Cannot yet use comm_style tiled with multi-mode comm");
		}

		/* ----------------------------------------------------------------------
		@@ -122,6 +136,7 @@ void CommTiled::setup()

		dimension = domain->dimension;
		int *periodicity = domain->periodicity;
		int ntypes = atom->ntypes;

		if (triclinic == 0) {
		prd = domain->prd;
		@@ -158,6 +173,18 @@ void CommTiled::setup()
		// set cutoff for comm forward and comm reverse
		// check that cutoff < any periodic box length

		if (mode == Comm::MULTI) {
		double *cuttype = neighbor->cuttype;
		double cut;
		for (i = 1; i <= ntypes; i++) {
		cut = 0.0;
		if (cutusermulti) cut = cutusermulti[i];
		cutghostmulti[i][0] = MAX(cut,cuttype[i]);
		cutghostmulti[i][1] = MAX(cut,cuttype[i]);
		cutghostmulti[i][2] = MAX(cut,cuttype[i]);
		}
		}

		double cut = get_comm_cutoff();
		if ((cut == 0.0) && (me == 0))
		error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
		@@ -173,6 +200,13 @@ void CommTiled::setup()
		cutghost[1] = cut * length1;
		length2 = h_inv[2];
		cutghost[2] = cut * length2;
		if (mode == Comm::MULTI){
		for (i = 1; i <= ntypes; i++) {
		cutghostmulti[i][0] *= length0;
		cutghostmulti[i][1] *= length1;
		cutghostmulti[i][2] *= length2;
		}
		}
		}

		if ((periodicity[0] && cutghost[0] > prd[0]) \|\|
		@@ -308,9 +342,11 @@ void CommTiled::setup()
		// = obox in other 2 dims
		// if sbox touches other proc's sub-box boundaries in lower dims,
		// extend sbox in those lower dims to include ghost atoms
		// single mode and multi mode

		double oboxlo[3],oboxhi[3],sbox[6];
		double oboxlo[3],oboxhi[3],sbox[6],sbox_multi[6];

		if (mode == Comm::SINGLE) {
		for (i = 0; i < noverlap; i++) {
		pbc_flag[iswap][i] = 0;
		pbc[iswap][i][0] = pbc[iswap][i][1] = pbc[iswap][i][2] =
		@@ -325,7 +361,6 @@ void CommTiled::setup()
		sbox[3] = MIN(oboxhi[0],hi1[0]);
		sbox[4] = MIN(oboxhi[1],hi1[1]);
		sbox[5] = MIN(oboxhi[2],hi1[2]);

		} else {
		pbc_flag[iswap][i] = 1;
		if (idir == 0) pbc[iswap][i][idim] = 1;
		@@ -334,7 +369,6 @@ void CommTiled::setup()
		if (idim == 1) pbc[iswap][i][5] = pbc[iswap][i][idim];
		if (idim == 2) pbc[iswap][i][4] = pbc[iswap][i][3] = pbc[iswap][i][idim];
		}

		sbox[0] = MAX(oboxlo[0],lo2[0]);
		sbox[1] = MAX(oboxlo[1],lo2[1]);
		sbox[2] = MAX(oboxlo[2],lo2[2]);
		@@ -347,7 +381,8 @@ void CommTiled::setup()
		sbox[idim] = sublo[idim];
		if (i < noverlap1)
		sbox[3+idim] = MIN(sbox[3+idim]+cutghost[idim],subhi[idim]);
		else sbox[3+idim] = MIN(sbox[3+idim]-prd[idim]+cutghost[idim],subhi[idim]);
		else
		sbox[3+idim] = MIN(sbox[3+idim]-prd[idim]+cutghost[idim],subhi[idim]);
		} else {
		if (i < noverlap1) sbox[idim] = MAX(sbox[idim]-cutghost[idim],sublo[idim]);
		else sbox[idim] = MAX(sbox[idim]+prd[idim]-cutghost[idim],sublo[idim]);
		@@ -365,6 +400,83 @@ void CommTiled::setup()

		memcpy(sendbox[iswap][i],sbox,6*sizeof(double));
		}
		}

		if (mode == Comm::MULTI) {
		for (i = 0; i < noverlap; i++) {
		pbc_flag[iswap][i] = 0;
		pbc[iswap][i][0] = pbc[iswap][i][1] = pbc[iswap][i][2] =
		pbc[iswap][i][3] = pbc[iswap][i][4] = pbc[iswap][i][5] = 0;

		(this->*box_other)(idim,idir,overlap[i],oboxlo,oboxhi);

		if (i < noverlap1) {
		sbox[0] = MAX(oboxlo[0],lo1[0]);
		sbox[1] = MAX(oboxlo[1],lo1[1]);
		sbox[2] = MAX(oboxlo[2],lo1[2]);
		sbox[3] = MIN(oboxhi[0],hi1[0]);
		sbox[4] = MIN(oboxhi[1],hi1[1]);
		sbox[5] = MIN(oboxhi[2],hi1[2]);
		} else {
		pbc_flag[iswap][i] = 1;
		if (idir == 0) pbc[iswap][i][idim] = 1;
		else pbc[iswap][i][idim] = -1;
		if (triclinic) {
		if (idim == 1) pbc[iswap][i][5] = pbc[iswap][i][idim];
		if (idim == 2) pbc[iswap][i][4] = pbc[iswap][i][3] = pbc[iswap][i][idim];
		}
		sbox[0] = MAX(oboxlo[0],lo2[0]);
		sbox[1] = MAX(oboxlo[1],lo2[1]);
		sbox[2] = MAX(oboxlo[2],lo2[2]);
		sbox[3] = MIN(oboxhi[0],hi2[0]);
		sbox[4] = MIN(oboxhi[1],hi2[1]);
		sbox[5] = MIN(oboxhi[2],hi2[2]);
		}

		for (int itype = 1; itype <= atom->ntypes; itype++) {
		sbox_multi[0] = sbox[0];
		sbox_multi[1] = sbox[1];
		sbox_multi[2] = sbox[2];
		sbox_multi[3] = sbox[3];
		sbox_multi[4] = sbox[4];
		sbox_multi[5] = sbox[5];
		if (idir == 0) {
		sbox_multi[idim] = sublo[idim];
		if (i < noverlap1)
		sbox_multi[3+idim] =
		MIN(sbox_multi[3+idim]+cutghostmulti[itype][idim],subhi[idim]);
		else
		sbox_multi[3+idim] =
		MIN(sbox_multi[3+idim]-prd[idim]+cutghostmulti[itype][idim],
		subhi[idim]);
		} else {
		if (i < noverlap1)
		sbox_multi[idim] =
		MAX(sbox_multi[idim]-cutghostmulti[itype][idim],sublo[idim]);
		else
		sbox_multi[idim] =
		MAX(sbox_multi[idim]+prd[idim]-cutghostmulti[itype][idim],
		sublo[idim]);
		sbox_multi[3+idim] = subhi[idim];
		}

		if (idim >= 1) {
		if (sbox_multi[0] == oboxlo[0])
		sbox_multi[0] -= cutghostmulti[itype][idim];
		if (sbox_multi[3] == oboxhi[0])
		sbox_multi[3] += cutghostmulti[itype][idim];
		}
		if (idim == 2) {
		if (sbox_multi[1] == oboxlo[1])
		sbox_multi[1] -= cutghostmulti[itype][idim];
		if (sbox_multi[4] == oboxhi[1])
		sbox_multi[4] += cutghostmulti[itype][idim];
		}

		memcpy(sendbox_multi[iswap][i][itype],sbox_multi,6*sizeof(double));
		}
		}
		}

		iswap++;
		}
		@@ -445,15 +557,15 @@ void CommTiled::setup()
		}
		}

		// reallocate MPI Requests and Statuses as needed
		// reallocate MPI Requests as needed

		int nmax = 0;
		for (i = 0; i < nswap; i++) nmax = MAX(nmax,nprocmax[i]);
		for (i = 0; i < dimension; i++) nmax = MAX(nmax,nexchprocmax[i]);
		if (nmax > maxreqstat) {
		maxreqstat = nmax;
		if (nmax > maxrequest) {
		maxrequest = nmax;
		delete [] requests;
		requests = new MPI_Request[maxreqstat];
		requests = new MPI_Request[maxrequest];
		}
		}

		@@ -815,7 +927,7 @@ void CommTiled::borders()
		// for x-dim swaps, check owned atoms
		// for yz-dim swaps, check owned and ghost atoms
		// store sent atom indices in sendlist for use in future timesteps
		// NOTE: assume SINGLE mode, add logic for MULTI mode later
		// single mode and multi mode

		x = atom->x;
		if (iswap % 2 == 0) nlast = atom->nlocal + atom->nghost;
		@@ -823,6 +935,8 @@ void CommTiled::borders()
		ncountall = 0;

		for (m = 0; m < nsendproc[iswap]; m++) {

		if (mode == Comm::SINGLE) {
		bbox = sendbox[iswap][m];
		xlo = bbox[0]; ylo = bbox[1]; zlo = bbox[2];
		xhi = bbox[3]; yhi = bbox[4]; zhi = bbox[5];
		@@ -861,7 +975,60 @@ void CommTiled::borders()
		sendnum[iswap][m] = ncount;
		smaxone = MAX(smaxone,ncount);
		ncountall += ncount;

		} else {

		int* type=atom->type;
		int itype;
		ncount = 0;

		if (!bordergroup) {
		for (i = 0; i < nlast; i++) {
		itype=type[i];
		bbox = sendbox_multi[iswap][m][itype];
		xlo = bbox[0]; ylo = bbox[1]; zlo = bbox[2];
		xhi = bbox[3]; yhi = bbox[4]; zhi = bbox[5];
		if (x[i][0] >= xlo && x[i][0] < xhi &&
		x[i][1] >= ylo && x[i][1] < yhi &&
		x[i][2] >= zlo && x[i][2] < zhi) {
		if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
		sendlist[iswap][m][ncount++] = i;
		}
		}
		} else {
		ngroup = atom->nfirst;
		for (i = 0; i < ngroup; i++) {
		itype=type[i];
		bbox = sendbox_multi[iswap][m][itype];
		xlo = bbox[0]; ylo = bbox[1]; zlo = bbox[2];
		xhi = bbox[3]; yhi = bbox[4]; zhi = bbox[5];
		if (x[i][0] >= xlo && x[i][0] < xhi &&
		x[i][1] >= ylo && x[i][1] < yhi &&
		x[i][2] >= zlo && x[i][2] < zhi) {
		if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
		sendlist[iswap][m][ncount++] = i;
		}
		}
		for (i = atom->nlocal; i < nlast; i++) {
		itype=type[i];
		bbox = sendbox_multi[iswap][m][itype];
		xlo = bbox[0]; ylo = bbox[1]; zlo = bbox[2];
		xhi = bbox[3]; yhi = bbox[4]; zhi = bbox[5];
		if (x[i][0] >= xlo && x[i][0] < xhi &&
		x[i][1] >= ylo && x[i][1] < yhi &&
		x[i][2] >= zlo && x[i][2] < zhi) {
		if (ncount == maxsendlist[iswap][m]) grow_list(iswap,m,ncount);
		sendlist[iswap][m][ncount++] = i;
		}
		}
		}

		sendnum[iswap][m] = ncount;
		smaxone = MAX(smaxone,ncount);
		ncountall += ncount;
		}
		}

		smaxall = MAX(smaxall,ncountall);

		// send sendnum counts to procs who recv from me except self
		@@ -914,9 +1081,8 @@ void CommTiled::borders()
		if (rmaxallsize_border > maxrecv) grow_recv(rmaxallsize_border);

		// swap atoms with other procs using pack_border(), unpack_border()
		// use Waitall() instead of Waitany() because calls to unpack_border()
		// must increment per-atom arrays in ascending order
		// For the same reason, sendself unpacks must occur after recvother unpacks
		// can use Waitany() because calls to unpack_border()
		// increment per-atom arrays as much as needed

		if (ghost_velocity) {
		if (recvother[iswap]) {
		@@ -932,13 +1098,6 @@ void CommTiled::borders()
		MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][m],0,world);
		}
		}
		if (recvother[iswap]) {
		MPI_Waitall(nrecv,requests,MPI_STATUS_IGNORE);
		for (m = 0; m < nrecv; m++)
		avec->unpack_border_vel(recvnum[iswap][m],firstrecv[iswap][m],
		&buf_recv[size_border*
		forward_recv_offset[iswap][m]]);
		}
		if (sendself[iswap]) {
		avec->pack_border_vel(sendnum[iswap][nsend],sendlist[iswap][nsend],
		buf_send,pbc_flag[iswap][nsend],
		@@ -946,6 +1105,14 @@ void CommTiled::borders()
		avec->unpack_border_vel(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
		buf_send);
		}
		if (recvother[iswap]) {
		for (i = 0; i < nrecv; i++) {
		MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE);
		avec->unpack_border_vel(recvnum[iswap][m],firstrecv[iswap][m],
		&buf_recv[size_border*
		forward_recv_offset[iswap][m]]);
		}
		}

		} else {
		if (recvother[iswap]) {
		@@ -961,19 +1128,20 @@ void CommTiled::borders()
		MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][m],0,world);
		}
		}
		if (recvother[iswap]) {
		MPI_Waitall(nrecv,requests,MPI_STATUS_IGNORE);
		for (m = 0; m < nrecv; m++)
		avec->unpack_border(recvnum[iswap][m],firstrecv[iswap][m],
		&buf_recv[size_border*
		forward_recv_offset[iswap][m]]);
		}
		if (sendself[iswap]) {
		avec->pack_border(sendnum[iswap][nsend],sendlist[iswap][nsend],
		buf_send,pbc_flag[iswap][nsend],pbc[iswap][nsend]);
		avec->unpack_border(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
		buf_send);
		}
		if (recvother[iswap]) {
		for (i = 0; i < nrecv; i++) {
		MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE);
		avec->unpack_border(recvnum[iswap][m],firstrecv[iswap][m],
		&buf_recv[size_border*
		forward_recv_offset[iswap][m]]);
		}
		}
		}

		// increment ghost atoms
		@@ -1476,11 +1644,9 @@ int CommTiled::exchange_variable(int n, double * /inbuf/, double & /outbuf*/

		void CommTiled::box_drop_brick(int idim, double lo, double hi, int &indexme)
		{
		// NOTE: this is not triclinic compatible
		// NOTE: these error messages are internal sanity checks
		// should not occur, can be removed at some point
		int dir;
		int index = -1;

		int index=-1,dir;
		if (hi[idim] == sublo[idim]) {
		index = myloc[idim] - 1;
		dir = -1;
		@@ -1913,6 +2079,7 @@ void CommTiled::allocate_swap(int n)
		pbc_flag = new int*[n];
		pbc = new int**[n];
		sendbox = new double**[n];
		sendbox_multi = new double***[n];
		maxsendlist = new int*[n];
		sendlist = new int**[n];

		@@ -1926,11 +2093,12 @@ void CommTiled::allocate_swap(int n)
		pbc_flag[i] = NULL;
		pbc[i] = NULL;
		sendbox[i] = NULL;
		sendbox_multi[i] = NULL;
		maxsendlist[i] = NULL;
		sendlist[i] = NULL;
		}

		maxreqstat = 0;
		maxrequest = 0;
		requests = NULL;

		for (int i = 0; i < n; i++) {
		@@ -1974,6 +2142,8 @@ void CommTiled::grow_swap_send(int i, int n, int nold)
		memory->create(pbc[i],n,6,"comm:pbc_flag");
		memory->destroy(sendbox[i]);
		memory->create(sendbox[i],n,6,"comm:sendbox");
		memory->destroy(sendbox_multi[i]);
		memory->create(sendbox_multi[i],n,atom->ntypes+1,6,"comm:sendbox_multi");

		delete [] maxsendlist[i];
		maxsendlist[i] = new int[n];
		@@ -2032,6 +2202,8 @@ void CommTiled::deallocate_swap(int n)
		delete [] pbc_flag[i];
		memory->destroy(pbc[i]);
		memory->destroy(sendbox[i]);
		memory->destroy(sendbox_multi[i]);

		delete [] maxsendlist[i];

		for (int j = 0; j < nprocmax[i]; j++) memory->destroy(sendlist[i][j]);
		@@ -2052,6 +2224,7 @@ void CommTiled::deallocate_swap(int n)
		delete [] pbc_flag;
		delete [] pbc;
		delete [] sendbox;
		delete [] sendbox_multi;
		delete [] maxsendlist;
		delete [] sendlist;

src/comm_tiled.h

+7 −3

Original line number	Diff line number	Diff line
		@@ -70,7 +70,6 @@ class CommTiled : public Comm {
		int **size_reverse_recv; // # of values to recv in each reverse swap/proc
		int **forward_recv_offset; // forward comm offsets in buf_recv per swap/proc
		int **reverse_recv_offset; // reverse comm offsets in buf_recv per swap/proc

		int ***sendlist; // list of atoms to send per swap/proc
		int **maxsendlist; // max size of send list per swap/proc
		int **pbc_flag; // general flag for sending atoms thru PBC
		@@ -78,6 +77,10 @@ class CommTiled : public Comm {

		double ***sendbox; // bounding box of atoms to send per swap/proc

		double **cutghostmulti; // cutghost on a per-type basis
		double ****sendbox_multi; // bounding box of atoms to send
		// per swap/proc for multi comm

		// exchange comm info, proc lists do not include self

		int *nexchproc; // # of procs to send/recv to/from in each dim
		@@ -92,7 +95,7 @@ class CommTiled : public Comm {
		int smaxall,rmaxall; // max size in atoms of any borders send/recv
		// for comm to all procs in one swap

		int maxreqstat; // max size of Request and Status vectors
		int maxrequest; // max size of Request vector
		MPI_Request *requests;

		struct RCBinfo {
		@@ -147,6 +150,7 @@ class CommTiled : public Comm {
		void grow_swap_send(int, int, int); // grow swap arrays for send and recv
		void grow_swap_recv(int, int);
		void deallocate_swap(int); // deallocate swap arrays

		};

		}

unittest/force-styles/test_pair_style.cpp

+5 −2

Original line number	Diff line number	Diff line
		@@ -1101,9 +1101,12 @@ TEST(PairStyle, single)

		// The single function in EAM is different from what we assume
		// here, therefore we have to skip testing those pair styles.
		if ((test_config.pair_style.substr(0, 3) == "eam") \|\|
		// Pair style colloid is also not compatible with this single tester
		if ((test_config.pair_style.substr(0, 7) == "colloid") \|\|
		(test_config.pair_style.substr(0, 3) == "eam") \|\|
		((test_config.pair_style.substr(0, 6) == "hybrid") &&
		(test_config.pair_style.find("eam") != std::string::npos))) {
		(test_config.pair_style.find("eam") != std::string::npos))
		) {
		if (!verbose) ::testing::internal::CaptureStdout();
		cleanup_lammps(lmp, test_config);
		if (!verbose) ::testing::internal::GetCapturedStdout();

Admin message