Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88751807
import_lammps.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Oct 20, 12:32
Size
21 KB
Mime Type
text/x-c++
Expires
Tue, Oct 22, 12:32 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21812455
Attached To
rLIBMULTISCALE LibMultiScale
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
Event Timeline
Log In to Comment