Homec4science
Diffusion LibMultiScale (features/components) Browse

libmultiscale/src/md/lammps/import_lammps.cc
b05d5f2aba78features/components

Last Change
Raw File

import_lammps.cc
View Options

/**
* @file import_lammps.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Wed Nov 20 07:29:13 2013
*
* @brief The programming hooks to let LM be aware of processor migrations
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#ifdef LIBMULTISCALE_USE_LAMMPS
#include "geometry.hh"
#include "domain_lammps.hh"
#include "reference_manager_lammps.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
using namespace LAMMPS_NS;
#include "memory.h"
#include "atom_vec.h"
#include "atom_vec_atomic.h"
#include "domain.h"
#include "fix.h"
#include "lammps_common.hh"
/* -------------------------------------------------------------------------- */
#include "import_lammps.hh"
/* -------------------------------------------------------------------------- */
enum Topology {
FACEXLEFT,
FACEXRIGHT,
FACEYLEFT,
FACEYRIGHT,
FACEZLEFT,
FACEZRIGHT,
EDGE1,
EDGE2,
EDGE3,
EDGE4,
EDGE5,
EDGE6,
EDGE7,
EDGE8,
EDGE9,
EDGE10,
EDGE11,
EDGE12,
CORNER1,
CORNER2,
CORNER3,
CORNER4,
CORNER5,
CORNER6,
CORNER7,
CORNER8
};
/* -------------------------------------------------------------------------- */
static const int n_neighs = 26;
static const int commPairs[n_neighs][2] = {
//faces pairs
{ FACEXLEFT ,FACEXRIGHT },
{ FACEXRIGHT,FACEXLEFT },
{ FACEYLEFT ,FACEYRIGHT },
{ FACEYRIGHT,FACEYLEFT },
{ FACEZLEFT ,FACEZRIGHT },
{ FACEZRIGHT,FACEZLEFT },
//edges pairs
{ EDGE1 , EDGE2 },
{ EDGE2 , EDGE1 },
{ EDGE3 , EDGE4 },
{ EDGE4 , EDGE3 },
{ EDGE5 , EDGE6 },
{ EDGE6 , EDGE5 },
{ EDGE7 , EDGE8 },
{ EDGE8 , EDGE7 },
{ EDGE9 , EDGE10 },
{ EDGE10, EDGE9 },
{ EDGE11, EDGE12 },
{ EDGE12, EDGE11 },
//corner pairs
{ CORNER1, CORNER2 },
{ CORNER2, CORNER1 },
{ CORNER3, CORNER4 },
{ CORNER4, CORNER3 },
{ CORNER5, CORNER6 },
{ CORNER6, CORNER5 },
{ CORNER7, CORNER8 },
{ CORNER8, CORNER7 }
};
/* -------------------------------------------------------------------------- */
int topologyCoords[n_neighs][3] = {
// faces
{-1, 0, 0}, //FACEXLEFT
{ 1, 0, 0}, //FACEXRIGHT
{ 0,-1, 0}, //FACEYLEFT
{ 0, 1, 0}, //FACEYRIGHT
{ 0, 0,-1}, //FACEZLEFT
{ 0, 0, 1}, //FACEZRIGHT
//edges
{ 1, 0,-1}, //EDGE1
{-1, 0,-1}, //EDGE2
{ 1, 0, 1}, //EDGE3
{-1, 0, 1}, //EDGE4
{-1,-1, 0}, //EDGE5
{ 1,-1, 0}, //EDGE6
{-1, 1, 0}, //EDGE7
{ 1, 1, 0}, //EDGE8
{ 0,-1,-1}, //EDGE9
{ 0, 1,-1}, //EDGE1
{ 0,-1, 1}, //EDGE1
{ 0, 1, 1}, //EDGE1
//corners
{-1,-1,-1},//CORNER1
{-1, 1,-1},//CORNER2
{-1, 1, 1},//CORNER3
{-1,-1, 1},//CORNER4
{1,-1,-1},//CORNER5
{1, 1,-1},//CORNER6
{1, 1, 1},//CORNER7
{1,-1, 1}//CORNER8
};
const char topologyStrings[n_neighs][255] = {
"FACEXLEFT",
"FACEXRIGHT",
"FACEYLEFT",
"FACEYRIGHT",
"FACEZLEFT",
"FACEZRIGHT",
"EDGE1",
"EDGE2",
"EDGE3",
"EDGE4",
"EDGE5",
"EDGE6",
"EDGE7",
"EDGE8",
"EDGE9",
"EDGE10",
"EDGE11",
"EDGE12",
"CORNER1",
"CORNER2",
"CORNER3",
"CORNER4",
"CORNER5",
"CORNER6",
"CORNER7",
"CORNER8"
};
/* -------------------------------------------------------------------------- */
#define stamp10 1000
#define stamp20 2000
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
void checkMPIComm(MPI_Status & status, int flag){
if (flag != MPI_SUCCESS || status.MPI_ERROR != MPI_SUCCESS) {
char errStr[MPI_MAX_ERROR_STRING];
int len;
MPI_Error_string(status.MPI_ERROR, errStr, &len);
LM_FATAL("MPI_ERROR = " << status.MPI_ERROR << "[" << errStr);
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::init(LAMMPS_NS::Domain & domain,
LAMMPS_NS::Atom & atom,
LAMMPS_NS::Comm & comm,
MPI_Comm & world,
int me, int nprocs){
old_nlocal = atom.nlocal;
VIEW_ATOM(RefLammps<Dim>);
this->atom = &atom;
this->domain = &domain;
this->comm = &comm;
this->world = world;
this->me = me;
this->nprocs = nprocs;
setPeriodic(domain.xperiodic,domain.yperiodic,domain.zperiodic);
setupTopology();
MPI_Barrier(world);
DUMP("before exchange nlocal = " << atom.nlocal,DBG_INFO);
if (domain.triclinic == 0) {
Xlo = domain.sublo;
Xhi = domain.subhi;
} else {
Xlo = domain.sublo_lamda;
Xhi = domain.subhi_lamda;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::registerSends(){
//clear the buffers
for (UInt i = 0; i < buf_send.size(); ++i)
buf_send[i].clear();
int nlocal = atom->nlocal;
int proc[3];
for (UInt i = 0; i < 3; ++i)
DUMP("box dimension " << i << " " << Xlo[i] << " " << Xhi[i],DBG_INFO);
DUMP(" nlocal = " << nlocal,DBG_INFO);
int at = 0;
DUMP("doing loop over atoms for migration",DBG_INFO);
while (at < nlocal) {
Real * X = atom->x[at];
#ifndef LM_OPTIMIZED
Real * X0 = atom->x0[at];
Real * V = atom->v[at];
Real * F = atom->f[at];
#endif
//setting actual coordinates
for (UInt k = 0; k < 3; ++k) {
proc[k] = 0;
// search the new processor
if (X[k] < Xlo[k]) proc[k] = -1;
if (X[k] >= Xhi[k]) proc[k] = 1;
//special treatment for pbc
if (periodic[k]) {
if (domain->triclinic == 0) {
//test if i am at the lower border
if (Xlo[k] == domain->boxlo[k])
if (X[k] >= Xhi[k] && (X[k] - Xhi[k]) > (Xhi[k]-Xlo[k]))
proc[k] = -1;
//test if i am at the upper border
if (Xhi[k] == domain->boxhi[k])
if (X[k] < Xlo[k] && (Xlo[k] - X[k]) > (Xhi[k]-Xlo[k]))
proc[k] = 1;
}
else {
//test if i am at the lower border
if (Xlo[k] == domain->boxlo_lamda[k])
if (X[k] >= Xhi[k] && (X[k] - Xhi[k]) > (Xhi[k]-Xlo[k]))
proc[k] = -1;
//test if i am at the upper border
if (Xhi[k] == domain->boxhi_lamda[k])
if (X[k] < Xlo[k] && (Xlo[k] - X[k]) > (Xhi[k]-Xlo[k]))
proc[k] = 1;
}
}
}
if (proc[0] == 0 && proc[1] == 0 && proc[2] == 0){
DUMP("No proc change for this atom " << at,DBG_ALL);
++at;
continue;
}
DUMP("migrating atom " << at << " nlocal << " << nlocal,DBG_ALL);
for (UInt k = 0; k < 3; ++k) {
DUMP("coord " << k << " "
<< X[k] << " " << X0[k] << " "
<< V[k] << " " << F[k] << " "
<< "[" << Xlo[k] << " - " << Xhi[k] << "] "
<< proc[k] << " " << proc[k],DBG_ALL);}
//DUMP("i: "<< at <<", X: " << X[0] << " is in proc " << proc[0],DBG_WARNING);
//DUMP("i: "<< at <<", Y: " << X[1] << " is in proc " << proc[1],DBG_WARNING);
//DUMP("i: "<< at <<", Z: " << X[2] << " is in proc " << proc[2],DBG_WARNING);
// build the sending buffers
for (int i = 0; i < n_neighs; ++i){
if (topologyCoords[i][0] == proc[0] &&
topologyCoords[i][1] == proc[1] &&
topologyCoords[i][2] == proc[2])
if (comProcs[i] != MPI_PROC_NULL) registerSend(at,comProcs[i]);
}
// fill the hole made by leaving atom
atom->avec->copy(nlocal-1,at);
// moveAtom(nlocal-1,at);
// decrease total number of atoms
nlocal--;
}
atom->nlocal = nlocal;
DUMP("nlocal is now " << atom->nlocal,DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::registerSend(UInt at, UInt toProc) {
#ifdef TRACE_ATOM
CHECK_TRACED(atom->x0[at]){
VIEW_ATOM(RefLammps<Dim>);
DUMP(ATOM_OUT(atom->x0[at])
<< " atom is leaving my proc for proc " << toProc,DBG_MESSAGE);
DUMP(ATOM_OUT(atom->x0[at])
<< " atom is leaving my proc for proc " << toProc,DBG_WARNING);
internal_index = -1;
}
#endif
sendAtom(at,toProc);
UInt nsend = buf_send[toProc].size();
buf_send[toProc].resize(nsend+packSize);
DUMP("buf_send size is now " << buf_send[toProc].size(),DBG_DETAIL);
UInt added = atom->avec->pack_exchange(at,&buf_send[toProc][nsend]);
if (added != packSize) LM_FATAL("internal error");
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::unPackAtoms(UInt procDe){
UInt m = 0;
std::vector<Real> & buf = buf_recv[procDe];
DUMP("unpacking atom " << buf_recv[procDe].size()/packSize << " atoms",DBG_INFO);
while (m < buf_recv[procDe].size()) {
m += atom->avec->unpack_exchange(&buf[m]);
acceptAtom(atom->nlocal-1,procDe);
{
UInt i = atom->nlocal-1;
Real** ptr = atom->x;
Real ** x0ptr = atom->x0;
Real ** vptr = atom->v;
SET_INTERNAL_TRACE_INDEX(x0ptr[i],i);
CHECK_TRACED(x0ptr[i])
{
DUMP("received traced atom from "
<< procDe,DBG_MESSAGE);
VIEW_ATOM(RefLammps<Dim>);
}
if ((ptr[i][0] < Xlo[0] || ptr[i][0] >= Xhi[0])
|| ptr[i][1] < Xlo[1] || ptr[i][1] >= Xhi[1]
|| ptr[i][2] < Xlo[2] || ptr[i][2] >= Xhi[2])
{
LM_FATAL("pb severe in migration process (received from "
<< procDe << " : atom i " << i
<< " x " << ptr[i][0] << " " << Xlo[0] << " " << Xhi[0] << " "
<< " y " << ptr[i][1] << " " << Xlo[1] << " " << Xhi[1] << " "
<< " z " << ptr[i][2] << " " << Xlo[2] << " " << Xhi[2] << " "
<< " big trouble an atom went far far away cannot "
<< "control anything now : its position is "
<< x0ptr[i][0] <<" "<< x0ptr[i][1] <<" "<< x0ptr[i][2]
<< " and its velocity is "
<< vptr[i][0] <<" "<< vptr[i][1] <<" "<< vptr[i][2]);
}
}
}
if (m != buf_recv[procDe].size()) LM_FATAL("not all atoms where unpacked");
buf_recv[procDe].clear();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::checkConsistency() {
MPI_Barrier(world);
int test=0;
MPI_Allreduce(&(atom->nlocal),&test,1,
MPI_INT, MPI_SUM, world);
DUMP("now nlocal = " << atom->nlocal
<< " nmax = " << atom->nmax,DBG_INFO);
if (test != atom->natoms){
DUMP("atoms were lost/created " << test << " != "
<< atom->natoms << " " << test-atom->natoms,DBG_MESSAGE);
DUMP("# local atoms is " << atom->nlocal,DBG_MESSAGE);
DUMP("old # local atoms is " << old_nlocal,DBG_MESSAGE);
//ref_manager.printBilan();
LM_FATAL("abort");
}
if (atom->nmax < atom->nlocal) LM_FATAL("pb with allocation of the atoms");
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::atomCommunications() {
//! pack the atoms needing a migration
registerSends();
//! exchange the atoms
atomCommunication();
//! unpack the data
UInt ncom = proc_comm.size();
for (UInt i = 0; i < ncom; ++i) unPackAtoms(proc_comm[i]);
//check constistency
checkConsistency();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void ImportLammps<Dim>::atomCommunication() {
std::map<int,MPI_Request> requests;
UInt ncom = proc_comm.size();
//! asynchronous send phase
for (UInt i = 0; i < ncom; ++i){
UInt toProc = proc_comm[i];
if (int(toProc) == me) continue;
int nbSend = buf_send[toProc].size();
DUMP("Send to proc " << toProc << " the " << nbSend << " Reals t:"
<< stamp20 + me,DBG_INFO);
// if (!nbSend) continue;
if (requests.count(toProc)) LM_FATAL("request to " << toProc << " already made");
MPI_Request & req = requests[toProc];
MPI_Isend(&buf_send[toProc][0], nbSend, MPI_DOUBLE, toProc, stamp20 + me,world, &req);
}
//! blocking receive phase
for (UInt i = 0; i < ncom; ++i){
UInt fromProc = proc_comm[i];
if (int(fromProc) == me) continue;
MPI_Status status;
status.MPI_ERROR = MPI_SUCCESS;
DUMP("probing from proc " << fromProc,DBG_INFO);
MPI_Probe(fromProc,stamp20 + fromProc,world,&status);
int nbRecv = 0;
MPI_Get_count(&status,MPI_DOUBLE,&nbRecv);
buf_recv[fromProc].resize(nbRecv);
DUMP("Receive from proc " << fromProc << " the " << nbRecv << " Reals t:"
<< stamp20 + fromProc,DBG_INFO);
// if (!nbRecv) continue;
int ret = MPI_Recv(&buf_recv[fromProc][0], nbRecv, MPI_DOUBLE, fromProc, stamp20 + fromProc,world, &status);
checkMPIComm(status,ret);
}
//! wait for sends to be complete
std::map<int,MPI_Request>::iterator it = requests.begin();
std::map<int,MPI_Request>::iterator end = requests.end();
while(it != end){
MPI_Status status;
status.MPI_ERROR = MPI_SUCCESS;
MPI_Request req = it->second;
int toProc = it->first;
DUMP("wait for com to " << toProc << " to be finished",DBG_INFO);
int ret = MPI_Wait(&req,&status);
if (buf_send[toProc].size()) {
checkMPIComm(status,ret);
}
DUMP("send to proc " << toProc << " finished ",DBG_INFO);
++it;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
int ImportLammps<Dim>::getProcNumber(int ind[3] , int offset[3]) {
/* Return the number of the processor located at the position i[0],i[1],i[2] */
/* It returns MPI_PROC_NULL if none */
/* It hanles periodic topology */
int index[3];
for (UInt i = 0; i < 3; ++i) index[i] = ind[i] + offset[i];
DUMP("ind " << ind[0] << " " << ind[1] << " " << ind[2],DBG_DETAIL);
DUMP("offset " << offset[0] << " " << offset[1] << " " << offset[2],DBG_DETAIL);
DUMP("periodic " << periodic[0] << " " << periodic[1] << " " << periodic[2],DBG_DETAIL);
DUMP("nbrProc " << nbrProc[0] << " " << nbrProc[1] << " " << nbrProc[2],DBG_DETAIL);
DUMP("index " << index[0] << " " << index[1] << " " << index[2],DBG_DETAIL);
//exclude out of the domain faces (negative side)
for (UInt i = 0; i < 3; ++i)
if (index[i] < 0 && !periodic[i]) return MPI_PROC_NULL;
else if (index[i] < 0) index[i] = nbrProc[i] - 1;
//exclude out of the domain faces (positive side)
for (UInt i = 0; i < 3; ++i)
if (index[i] >= nbrProc[i] && !periodic[i]) return MPI_PROC_NULL;
else if (index[i] >= nbrProc[i]) index[i] = 0;
DUMP("index " << index[0] << " " << index[1] << " " << index[2],DBG_DETAIL);
return proc_matrix[index];
}
/* -------------------------------------------------------------------------- */
//!setting processor topology with a stamp style
template <UInt Dim>
void ImportLammps<Dim>::setupTopology() {
int reorder = 0;
// get number of processors in x,y and z directions
for (UInt i = 0; i < 3; ++i) nbrProc[i] = comm->procgrid[i];
// MPI communicator in cartesian coordinates
MPI_Comm cartesian;
// creation of the topology
MPI_Cart_create(world, 3, nbrProc, periodic, reorder, &cartesian);
// fill the processor repartition matrix
proc_matrix.resize(nbrProc);
int coords[3];
for (int n = 0; n < nprocs; n++) {
MPI_Cart_coords(cartesian, n, 3, coords);
proc_matrix[coords] = n;
DUMP("proc " << n << " has cartesian coordinates : "
<< coords[0] << " "
<< coords[1] << " "
<< coords[2],DBG_INFO);
}
/* coordonnes du processeur */
MPI_Cart_coords(cartesian, me, 3, coords);
DUMP("my coordinates are "
<< coords[0] << " "
<< coords[1] << " "
<< coords[2],DBG_INFO);
MPI_Comm_free(&cartesian);
// initialize comprocs
for (int i = 0; i < n_neighs; ++i){
DUMP("get " << topologyStrings[i] << " proc number",DBG_DETAIL);
comProcs[i] = getProcNumber(coords,topologyCoords[i]);
DUMP(topologyStrings[i] << " " << comProcs[i],DBG_INFO);
}
for (UInt i = 0; i < 3; ++i)
DUMP( "nbrProc[" << i << "]" << nbrProc[i],DBG_INFO);
//initialize com buffers array
buf_send.resize(nprocs);
buf_recv.resize(nprocs);
// init list of processors to discuss with
std::set<int> proc_comm_with;
for (int i = 0; i < n_neighs; ++i) proc_comm_with.insert(comProcs[i]);
proc_comm.clear();
std::set<int>::iterator it = proc_comm_with.begin();
std::set<int>::iterator end = proc_comm_with.end();
while(it != end){
if (*it != MPI_PROC_NULL) {
proc_comm.push_back(*it);
DUMP("I should communicate with proc " << *it,DBG_DETAIL);
}
++it;
}
DUMP("I should communicate with " << proc_comm.size() << " processors",DBG_DETAIL);
// compute the packSize
double tmp[1000];
packSize = atom->avec->pack_exchange(0,tmp);
DUMP("packSize is " << packSize,DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template class ImportLammps<2>;
template class ImportLammps<3>;
ImportLammpsInterface * ImportLammpsInterface::import = NULL;
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
using namespace libmultiscale;
/* ---------------------------------------------------------------------- */
void AtomVecAtomic::copy(int i, int j)
{
tag[j] = tag[i];
type[j] = type[i];
mask[j] = mask[i];
image[j] = image[i];
x[j][0] = x[i][0];
x[j][1] = x[i][1];
x[j][2] = x[i][2];
x0[j][0] = x0[i][0];
x0[j][1] = x0[i][1];
x0[j][2] = x0[i][2];
v[j][0] = v[i][0];
v[j][1] = v[i][1];
v[j][2] = v[i][2];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
ImportLammpsInterface::getImport().moveAtom(i,j);
}
/* -------------------------------------------------------------------------- */
void Comm::exchange() {
MPI_Barrier(world);
DUMP("migration coms",DBG_INFO);
ImportLammpsInterface::getImport().init(*domain,*atom,*this,world,this->me,this->nprocs);
// clear global->local map since atoms move & new ghosts are created
if (map_style) atom->map_clear();
DUMP("after construction of sending buffers nlocal = " << atom->nlocal
<< " nmax = " << atom->nmax,DBG_INFO);
// send/receive communication to other procs
ImportLammpsInterface::getImport().atomCommunications();
MPI_Barrier(world);
ImportLammpsInterface::getImport().updateRefSubSets();
MPI_Barrier(world);
DUMP("migration all done for me !",DBG_INFO);
}
/* -------------------------------------------------------------------------- */
// cleaned
/* -------------------------------------------------------------------------- */
void AtomVecAtomic::create_atom(int itype,Real *coord) {
// add the atom to my list of atoms
if (!ImportLammpsInterface::getImport().getGeom(itype).
contains(coord[0],coord[1],coord[2]))
return;
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = 0;
type[nlocal] = itype;
for (UInt i = 0; i < 3; ++i) {
x[nlocal][i] = coord[i];
x0[nlocal][i] = coord[i];
v[nlocal][i] = 0.0;
}
mask[nlocal] = 1;
image[nlocal] = (512 << 20) | (512 << 10) | 512;
Real * X = x[nlocal];
Real * xlo = domain->sublo;
Real * xhi = domain->subhi;
for (UInt i = 0; i < 3; ++i) {
if (X[i] < xlo[i]){
LM_FATAL("cannot accept to create initially an atom outside"
<< " of the initial box" <<
"X[" << i << "] = " << X[i] << " < xlo = " << xlo[i]);
}
if (X[i] >= xhi[i]){
LM_FATAL("cannot accept to create initially an atom outside"
<< " of the initial box" <<
"X[" << i << "] = " << X[i] << " >= xhi = "<< xhi[i]);
}
}
atom->nlocal++;
}
/* -------------------------------------------------------------------------- */
#endif
c4science · Help