Newer
Older
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, haktulga@cs.purdue.edu
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
Kurt A. O'Hearn
committed
#include "reax_types.h"
#include "grid.h"
#include "reset_tools.h"
#include "tool_box.h"
#include "vector.h"
Kurt A. O'Hearn
committed
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#if defined(NEUTRAL_TERRITORY)
void Setup_NT_Comm( reax_system * const system, control_params * const control,
mpi_datatypes * const mpi_data )
{
int i, d;
real bndry_cut;
neighbor_proc *nbr_pr;
simulation_box *my_box;
ivec nbr_coords, nbr_recv_coords;
ivec r[12] = {
{0, 0, -1}, // -z
{0, 0, +1}, // +z
{0, -1, 0}, // -y
{-1, -1, 0}, // -x-y
{-1, 0, 0}, // -x
{-1, +1, 0}, // -x+y
{0, 0, +1}, // +z
{0, 0, -1}, // -z
{0, +1, 0}, // +y
{+1, +1, 0}, // +x+y
{+1, 0, 0}, // +x
{+1, -1, 0} // +x-y
};
my_box = &system->my_box;
bndry_cut = system->bndry_cuts.ghost_cutoff;
system->num_nt_nbrs = MAX_NT_NBRS;
/* identify my neighbors */
for ( i = 0; i < system->num_nt_nbrs; ++i )
{
nbr_pr = &system->my_nt_nbrs[i];
ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &nbr_pr->rank );
/* set the rank of the neighbor processor in the receiving direction */
ivec_Sum( nbr_recv_coords, system->my_coords, r[i + 6] ); /* actual nbr coords */
MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_recv_coords, &nbr_pr->receive_rank );
for ( d = 0; d < 3; ++d )
{
/* determine the boundary area with this nbr */
if ( r[i][d] < 0 )
{
nbr_pr->bndry_min[d] = my_box->min[d];
nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
}
else if ( r[i][d] > 0 )
{
nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
nbr_pr->bndry_max[d] = my_box->max[d];
}
else
{
nbr_pr->bndry_min[d] = my_box->min[d];
nbr_pr->bndry_max[d] = my_box->max[d];
}
/* determine if it is a periodic neighbor */
if ( nbr_coords[d] < 0 )
{
nbr_pr->prdc[d] = -1;
}
else if ( nbr_coords[d] >= control->procs_by_dim[d] )
{
nbr_pr->prdc[d] = 1;
}
else
{
nbr_pr->prdc[d] = 0;
}
}
}
}
static int Sort_Neutral_Territory( reax_system *system, int dir, mpi_out_data *out_bufs, int write )
{
int i, cnt;
reax_atom *atoms;
neighbor_proc *nbr_pr;
cnt = 0;
atoms = system->my_atoms;
/* place each atom into the appropriate outgoing list */
nbr_pr = &( system->my_nt_nbrs[dir] );
for ( i = 0; i < system->n; ++i )
{
if ( nbr_pr->bndry_min[0] <= atoms[i].x[0]
&& atoms[i].x[0] < nbr_pr->bndry_max[0]
&& nbr_pr->bndry_min[1] <= atoms[i].x[1]
&& atoms[i].x[1] < nbr_pr->bndry_max[1]
&& nbr_pr->bndry_min[2] <= atoms[i].x[2]
&& atoms[i].x[2] < nbr_pr->bndry_max[2] )
{
if ( write )
{
out_bufs[dir].index[out_bufs[dir].cnt] = i;
out_bufs[dir].cnt++;
}
else
{
cnt++;
}
}
}
return cnt;
}
static void Init_Neutral_Territory( reax_system* system, mpi_datatypes *mpi_data )
{
int d, end, cnt;
mpi_out_data *out_bufs;
MPI_Comm comm;
MPI_Request req;
MPI_Status stat;
neighbor_proc *nbr;
Reset_Out_Buffers( mpi_data->out_nt_buffers, system->num_nt_nbrs );
comm = mpi_data->comm_mesh3D;
out_bufs = mpi_data->out_nt_buffers;
cnt = 0;
end = system->n;
for ( d = 0; d < 6; ++d )
{
nbr = &system->my_nt_nbrs[d];
Sort_Neutral_Territory( system, d, out_bufs, 1 );
MPI_Irecv( &cnt, 1, MPI_INT, nbr->receive_rank, d, comm, &req );
MPI_Send( &out_bufs[d].cnt, 1, MPI_INT, nbr->rank, d, comm );
MPI_Wait( &req, &stat );
if ( mpi_data->in_nt_buffer[d] == NULL )
{
nbr->est_recv = MAX( SAFER_ZONE_NT * cnt, MIN_SEND );
mpi_data->in_nt_buffer[d] = smalloc( nbr->est_recv * sizeof(real),
"Init_Neural_Territory::mpi_data->in_nt_buffer[d]", comm );
}
nbr = &system->my_nt_nbrs[d];
nbr->atoms_str = end;
nbr->atoms_cnt = cnt;
end += cnt;
}
}
void Estimate_NT_Atoms( reax_system * const system, mpi_datatypes * const mpi_data )
{
int d;
mpi_out_data *out_bufs;
neighbor_proc *nbr;
out_bufs = mpi_data->out_nt_buffers;
for ( d = 0; d < 6; ++d )
{
/* count the number of atoms in each processor's outgoing list */
nbr = &system->my_nt_nbrs[d];
nbr->est_send = Sort_Neutral_Territory( system, d, out_bufs, 0 );
/* estimate the space needed based on the count above */
nbr->est_send = MAX( MIN_SEND, nbr->est_send * SAFER_ZONE_NT );
/* allocate the estimated space */
out_bufs[d].index = scalloc( 2 * nbr->est_send, sizeof(int),
Kurt A. O'Hearn
committed
"Estimate_NT_Atoms::out_bufs[d].index", MPI_COMM_WORLD );
out_bufs[d].out_atoms = scalloc( 2 * nbr->est_send, sizeof(real),
Kurt A. O'Hearn
committed
"Estimate_NT_Atoms::out_bufs[d].out_atoms", MPI_COMM_WORLD );
/* sort the atoms to their outgoing buffers */
// TODO: to call or not to call?
//Sort_Neutral_Territory( system, d, out_bufs, 1 );
}
}
#endif
Kurt A. O'Hearn
committed
void Check_MPI_Error( int code, const char * const msg )
{
char err_msg[MPI_MAX_ERROR_STRING];
int len;
if ( code != MPI_SUCCESS )
{
MPI_Error_string( code, err_msg, &len );
Kurt A. O'Hearn
committed
fprintf( stderr, "[ERROR] MPI error code %d, from %s\n",
Kurt A. O'Hearn
committed
fprintf( stderr, " [INFO] MPI error message: %s\n", err_msg );
MPI_Abort( MPI_COMM_WORLD, RUNTIME_ERROR );
Kurt A. O'Hearn
committed
void Setup_Comm( reax_system * const system, control_params * const control,
mpi_datatypes * const mpi_data )
Kurt A. O'Hearn
committed
const real bndry_cut = system->bndry_cuts.ghost_cutoff;
Kurt A. O'Hearn
committed
simulation_box * const my_box = &system->my_box;
{ -1, 0, 0}, { 1, 0, 0}, // -x, +x
{0, -1, 0}, {0, 1, 0}, // -y, +y
{0, 0, -1}, {0, 0, 1}, // -z, +z
/* identify my neighbors */
system->num_nbrs = MAX_NBRS;
for ( i = 0; i < system->num_nbrs; ++i )
{
ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
nbr_pr = &system->my_nbrs[i];
Kurt A. O'Hearn
committed
MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &nbr_pr->rank );
for ( d = 0; d < 3; ++d )
{
/* determine the boundary area with this nbr */
if ( r[i][d] < 0 )
{
nbr_pr->bndry_min[d] = my_box->min[d];
nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
}
else if ( r[i][d] > 0 )
{
nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
nbr_pr->bndry_max[d] = my_box->max[d];
}
else
{
nbr_pr->bndry_min[d] = nbr_pr->bndry_max[d] = NEG_INF;
}
/* determine if it is a periodic neighbor */
if ( nbr_coords[d] < 0 )
else if ( nbr_coords[d] >= control->procs_by_dim[d] )
fprintf( stderr, "p%d-nbr%d: r[%2d %2d %2d] -> c[%2d %2d %2d] -> rank=%d\n",
system->my_rank, i, r[i][0], r[i][1], r[i][2],
nbr_coords[0], nbr_coords[1], nbr_coords[2], nbr_pr->rank );
Kurt A. O'Hearn
committed
void Update_Comm( reax_system * const system )
Kurt A. O'Hearn
committed
const real bndry_cut = system->bndry_cuts.ghost_cutoff;
Kurt A. O'Hearn
committed
simulation_box * const my_box = &system->my_box;
{ -1, 0, 0}, { 1, 0, 0}, // -x, +x
{0, -1, 0}, {0, 1, 0}, // -y, +y
{0, 0, -1}, {0, 0, 1}, // -z, +z
/* identify my neighbors */
for ( i = 0; i < system->num_nbrs; ++i )
{
nbr_pr = &system->my_nbrs[i];
/* determine the boundary area with this nbr */
if ( r[i][d] < 0 )
{
nbr_pr->bndry_min[d] = my_box->min[d];
nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
}
else if ( r[i][d] > 0 )
{
nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
nbr_pr->bndry_max[d] = my_box->max[d];
}
else
{
nbr_pr->bndry_min[d] = nbr_pr->bndry_max[d] = NEG_INF;
}
}
/********************* ATOM TRANSFER ***********************/
/***************** PACK & UNPACK ATOMS *********************/
Kurt A. O'Hearn
committed
static void Pack_MPI_Atom( mpi_atom * const matm, const reax_atom * const ratm, int i )
matm->orig_id = ratm->orig_id;
matm->type = ratm->type;
matm->num_bonds = ratm->num_bonds;
matm->num_hbonds = ratm->num_hbonds;
strncpy( matm->name, ratm->name, MAX_ATOM_NAME_LEN );
rvec_Copy( matm->x, ratm->x );
rvec_Copy( matm->v, ratm->v );
rvec_Copy( matm->f_old, ratm->f_old );
memcpy( matm->s, ratm->s, sizeof(rvec4) ); //rvec_Copy( matm->s, ratm->s );
memcpy( matm->t, ratm->t, sizeof(rvec4) ); //rvec_Copy( matm->t, ratm->t );
Kurt A. O'Hearn
committed
static void Unpack_MPI_Atom( reax_atom * const ratm, const mpi_atom * const matm )
ratm->orig_id = matm->orig_id;
ratm->type = matm->type;
ratm->num_bonds = matm->num_bonds;
ratm->num_hbonds = matm->num_hbonds;
strncpy( ratm->name, matm->name, MAX_ATOM_NAME_LEN );
rvec_Copy( ratm->x, matm->x );
rvec_Copy( ratm->v, matm->v );
rvec_Copy( ratm->f_old, matm->f_old );
memcpy( ratm->s, matm->s, sizeof(rvec4) ); //rvec_Copy( ratm->s, matm->s );
memcpy( ratm->t, matm->t, sizeof(rvec4) ); //rvec_Copy( ratm->t, matm->t );
}
/*********************** SORTER **************************/
Kurt A. O'Hearn
committed
static void Sort_Transfer_Atoms( reax_system * const system, int start, int end,
int dim, mpi_out_data * const out_bufs )
Kurt A. O'Hearn
committed
reax_atom * const atoms = system->my_atoms;
simulation_box * const my_box = &system->my_box;
Kurt A. O'Hearn
committed
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d sort_transfers: start=%d end=%d dim=%d starting...\n",
Kurt A. O'Hearn
committed
system->my_rank, start, end, dim );
/* place each atom into the appropriate outgoing list */
for ( i = start; i < end; ++i )
{
for ( d = dim; d < 3; ++d )
if ( atoms[i].x[d] < my_box->min[d] )
{
out_cnt = out_bufs[2 * d].cnt++;
Kurt A. O'Hearn
committed
out_buf = out_bufs[2 * d].out_atoms;
Pack_MPI_Atom( out_buf + out_cnt, atoms + i, i );
atoms[i].orig_id = -1;
break;
}
else if ( atoms[i].x[d] >= my_box->max[d] )
{
out_cnt = out_bufs[2 * d + 1].cnt++;
Kurt A. O'Hearn
committed
out_buf = out_bufs[2 * d + 1].out_atoms;
Pack_MPI_Atom( out_buf + out_cnt, atoms + i, i );
atoms[i].orig_id = -1;
break;
}
if ( out_bufs[d].cnt )
{
fprintf( stderr, "p%d to p%d(nbr%d) # of transfers = %d\n",
system->my_rank, system->my_nbrs[d].rank, d, out_bufs[d].cnt );
Kurt A. O'Hearn
committed
Kurt A. O'Hearn
committed
Kurt A. O'Hearn
committed
{
fprintf( stderr, "p%d to p%d: transfer atom%d [%.3f %.3f %.3f]\n",
system->my_rank, system->my_nbrs[d].rank, out_buf[i].imprt_id,
out_buf[i].x[0], out_buf[i].x[1], out_buf[i].x[2] );
Kurt A. O'Hearn
committed
}
}
//fprintf( stderr, "p%d sort_transfers: start=%d end=%d dim=%d done!\n",
// system->my_rank, start, end, dim );
}
/*********************** UNPACKER **************************/
Kurt A. O'Hearn
committed
static void Unpack_Transfer_Message( reax_system * const system, int end, void * const dummy,
int cnt, neighbor_proc * const nbr, int dim )
Kurt A. O'Hearn
committed
reax_atom * const dest = system->my_atoms + end;
Kurt A. O'Hearn
committed
mpi_atom * const src = dummy;
/* adjust coordinates of recved atoms if nbr is a periodic one */
if ( nbr->prdc[dim] )
{
dx = nbr->prdc[dim] * system->big_box.box_norms[dim];
Kurt A. O'Hearn
committed
}
/************** EXCHANGE BOUNDARY ATOMS *****************/
/************ PACK & UNPACK BOUNDARY ATOMS **************/
Kurt A. O'Hearn
committed
static void Pack_Boundary_Atom( boundary_atom * const matm, const reax_atom * const ratm, int i )
matm->orig_id = ratm->orig_id;
matm->type = ratm->type;
matm->num_bonds = ratm->num_bonds;
matm->num_hbonds = ratm->num_hbonds;
rvec_Copy( matm->x, ratm->x );
Kurt A. O'Hearn
committed
static void Unpack_Boundary_Atom( reax_atom * const ratm, const boundary_atom * const matm )
ratm->orig_id = matm->orig_id;
ratm->imprt_id = matm->imprt_id;
ratm->type = matm->type;
ratm->num_bonds = matm->num_bonds;
ratm->num_hbonds = matm->num_hbonds;
//ratm->renumber = offset + matm->imprt_id;
rvec_Copy( ratm->x, matm->x );
}
/*********************** SORTER **************************/
Kurt A. O'Hearn
committed
static void Sort_Boundary_Atoms( reax_system * const system, int start, int end,
int dim, mpi_out_data * const out_bufs )
Kurt A. O'Hearn
committed
const reax_atom * const atoms = system->my_atoms;
fprintf( stderr, "p%d sort_exchange: start=%d end=%d dim=%d starting...\n",
Kurt A. O'Hearn
committed
system->my_rank, start, end, dim );
/* place each atom into the appropriate outgoing list */
for ( i = start; i < end; ++i )
nbr_pr = &system->my_nbrs[p];
if ( nbr_pr->bndry_min[d] <= atoms[i].x[d] &&
atoms[i].x[d] < nbr_pr->bndry_max[d] )
{
out_cnt = out_bufs[p].cnt++;
out_bufs[p].index[out_cnt] = i;
Kurt A. O'Hearn
committed
out_buf = out_bufs[p].out_atoms;
Pack_Boundary_Atom( out_buf + out_cnt, atoms + i, i );
}
}
fprintf( stderr, "p%d to p%d(nbr%d) # of exchanges to send = %d\n",
system->my_rank, system->my_nbrs[i].rank, i,
out_bufs[i].cnt );
fprintf( stderr, "p%d sort_exchange: start=%d end=%d dim=%d done!\n",
system->my_rank, start, end, dim );
#endif
Kurt A. O'Hearn
committed
void Estimate_Boundary_Atoms( reax_system * const system, int start, int end,
int d, mpi_out_data * const out_bufs )
Kurt A. O'Hearn
committed
const reax_atom * const atoms = system->my_atoms;
Kurt A. O'Hearn
committed
neighbor_proc * const nbr1 = &system->my_nbrs[2 * d];
neighbor_proc * const nbr2 = &system->my_nbrs[2 * d + 1];
neighbor_proc *nbr_pr;
fprintf( stderr, "p%d estimate_exchange: start=%d end=%d dim=%d starting.\n",
Kurt A. O'Hearn
committed
system->my_rank, start, end, d );
nbr1->est_send = 0;
nbr2->est_send = 0;
/* count the number of atoms in each processor's outgoing list */
for ( i = 0; i < end; ++i )
{
if ( nbr1->bndry_min[d] <= atoms[i].x[d] && atoms[i].x[d] < nbr1->bndry_max[d] )
if ( nbr2->bndry_min[d] <= atoms[i].x[d] && atoms[i].x[d] < nbr2->bndry_max[d] )
}
/* estimate the space based on the count above */
nbr1->est_send = MAX( MIN_SEND, nbr1->est_send * SAFER_ZONE );
nbr2->est_send = MAX( MIN_SEND, nbr2->est_send * SAFER_ZONE );
fprintf( stderr, "p%d estimate_exchange: end=%d dim=%d est1=%d est2=%d\n",
Kurt A. O'Hearn
committed
system->my_rank, end, d, nbr1->est_send, nbr2->est_send );
/* allocate the estimated space */
for ( p = 2 * d; p < 2 * d + 2; ++p )
{
Kurt A. O'Hearn
committed
nbr_pr = &system->my_nbrs[p];
out_bufs[p].index = scalloc( 2 * nbr_pr->est_send, sizeof(int),
Kurt A. O'Hearn
committed
"Estimate_Boundary_Atoms::mpibuf:index" );
out_bufs[p].out_atoms = scalloc( 2 * nbr_pr->est_send, sizeof(boundary_atom),
Kurt A. O'Hearn
committed
"Estimate_Boundary_Atoms::mpibuf:out_atoms" );
/* sort the atoms to their outgoing buffers */
for ( i = 0; i < end; ++i )
/* check if atom is outbound to another processor
* in either direction of the dimension under consideration */
Kurt A. O'Hearn
committed
nbr_pr = &system->my_nbrs[p];
if ( nbr_pr->bndry_min[d] <= atoms[i].x[d] &&
atoms[i].x[d] < nbr_pr->bndry_max[d] )
{
out_cnt = out_bufs[p].cnt;
Kurt A. O'Hearn
committed
out_buf = out_bufs[p].out_atoms;
Pack_Boundary_Atom( out_buf + out_cnt, atoms + i, i );
for ( p = 2 * d; p < 2 * d + 2; ++p )
{
for ( i = 0; i < out_bufs[p].cnt; ++i )
{
fprintf( stderr, "p%d: out_bufs[%d].index[%d] = %d\n",
system->my_rank, p, i, out_bufs[p].index[i] );
fprintf( stderr, " p%d: atom %6d, x[0] = %10.4f, x[1] = %10.4f, x[2] = %10.4f\n",
system->my_rank,
((boundary_atom *)(out_bufs[p].out_atoms))[i].orig_id,
((boundary_atom *)(out_bufs[p].out_atoms))[i].x[0],
((boundary_atom *)(out_bufs[p].out_atoms))[i].x[1],
((boundary_atom *)(out_bufs[p].out_atoms))[i].x[2] );
}
}
fprintf( stderr, "p%d estimate_exchange: end=%d dim=%d done!\n",
system->my_rank, end, d );
Kurt A. O'Hearn
committed
static void Estimate_Init_Storage( int me, neighbor_proc * const nbr1, neighbor_proc * const nbr2,
int d, int * const max, int * const nrecv, mpi_datatypes * const mpi_data, MPI_Comm comm )
/* first exchange the estimates, then allocate buffers */
ret = MPI_Irecv( &nbr1->est_recv, 1, MPI_INT, nbr1->rank, 2 * d + 1, comm, &req1 );
Check_MPI_Error( ret, "Estimate_Init_Storage::MPI_Irecv::nbr1" );
ret = MPI_Irecv( &nbr2->est_recv, 1, MPI_INT, nbr2->rank, 2 * d, comm, &req2 );
Check_MPI_Error( ret, "Estimate_Init_Storage::MPI_Irecv::nbr2" );
ret = MPI_Send( &nbr1->est_send, 1, MPI_INT, nbr1->rank, 2 * d, comm );
Check_MPI_Error( ret, "Estimate_Init_Storage::MPI_Send::nbr1" );
ret = MPI_Send( &nbr2->est_send, 1, MPI_INT, nbr2->rank, 2 * d + 1, comm );
Check_MPI_Error( ret, "Estimate_Init_Storage::MPI_Send::nbr2" );
ret = MPI_Wait( &req1, &stat1 );
Check_MPI_Error( ret, "Estimate_Init_Storage::MPI_Wait::nbr1" );
ret = MPI_Wait( &req2, &stat2 );
Check_MPI_Error( ret, "Estimate_Init_Storage::MPI_Wait::nbr2" );
nrecv[2 * d] = nbr1->est_recv;
nrecv[2 * d + 1] = nbr2->est_recv;
new_max = MAX( nbr1->est_recv, nbr2->est_recv );
fprintf( stderr, "p%d-p%d(nbr%d) est_send=%d est_recv=%d\n",
me, nbr1->rank, 2 * d, nbr1->est_send, nbr1->est_recv );
fprintf( stderr, "p%d-p%d(nbr%d) est_send=%d est_recv=%d\n",
me, nbr2->rank, 2 * d + 1, nbr2->est_send, nbr2->est_recv );
fprintf( stderr, "max=%d new_max=%d\n", *max, new_max );
Kurt A. O'Hearn
committed
if ( mpi_data->in1_buffer != NULL )
Kurt A. O'Hearn
committed
sfree( mpi_data->in1_buffer, "Estimate_Init_Storage::mpi_data->in1_buffer" );
Kurt A. O'Hearn
committed
if ( mpi_data->in2_buffer != NULL )
Kurt A. O'Hearn
committed
sfree( mpi_data->in2_buffer, "Estimate_Init_Storage::mpi_data->in2_buffer" );
Kurt A. O'Hearn
committed
mpi_data->in1_buffer = smalloc( 2 * new_max * MAX3( sizeof(mpi_atom), sizeof(boundary_atom), sizeof(rvec) ),
Kurt A. O'Hearn
committed
"Estimate_Init_Storage::mpi_data->in1_buffer" );
mpi_data->in2_buffer = smalloc( 2 * new_max * MAX3( sizeof(mpi_atom), sizeof(boundary_atom), sizeof(rvec) ),
Kurt A. O'Hearn
committed
"Estimate_Init_Storage::mpi_data->in2_buffer" );
}
/*********************** UNPACKER **************************/
Kurt A. O'Hearn
committed
static void Unpack_Exchange_Message( reax_system * const system, int end, void * const dummy,
int cnt, neighbor_proc * const nbr, int dim )
Kurt A. O'Hearn
committed
const boundary_atom * const src = dummy;
Kurt A. O'Hearn
committed
reax_atom * const dest = &system->my_atoms[end];
Kurt A. O'Hearn
committed
Unpack_Boundary_Atom( &dest[i], &src[i] );
#if defined(DEBUG_FOCUS)
for ( i = end; i < end + cnt; ++i )
{
fprintf( stderr, "UNPACK p%d: d = %d, atom %d, x[0] = %10.4f, x[1] = %10.4f, x[2] = %10.4f\n",
system->my_rank, dim, i,
system->my_atoms[i].x[0],
system->my_atoms[i].x[1],
system->my_atoms[i].x[2] );
}
#endif
/* record the atoms recv'd from this nbr */
nbr->atoms_str = end;
nbr->atoms_cnt = cnt;
nbr->est_recv = MAX( (int)(cnt * SAFER_ZONE), MIN_SEND );
/* update max_recv to make sure that we reallocate at the right time */
if ( cnt > system->max_recved )
/* adjust coordinates of recved atoms if nbr is a periodic one */
if ( nbr->prdc[dim] )
{
dx = nbr->prdc[dim] * system->big_box.box_norms[dim];
#if defined(DEBUG_FOCUS)
fprintf( stderr, "UNPACK p%d: dim = %d, dx = %f\n",
system->my_rank, dim, dx );
#endif
Kurt A. O'Hearn
committed
void Unpack_Estimate_Message( reax_system * const system, int end, void * const dummy,
int cnt, neighbor_proc * const nbr, int dim )
fprintf( stderr, "p%d-p%d unpack_estimate: end=%d cnt=%d - unpacking\n",
system->my_rank, nbr->rank, end, cnt );
system->my_atoms = srealloc( system->my_atoms, sizeof(reax_atom) * (end + cnt),
"Unpack_Estimate_Message::system:my_atoms" );
Unpack_Exchange_Message( system, end, dummy, cnt, nbr, dim );
fprintf( stderr, "p%d-p%d unpack_estimate: end=%d cnt=%d - done\n",
system->my_rank, nbr->rank, end, cnt );
#endif
}
/**************** UPDATE ATOM POSITIONS *******************/
/**************** PACK POSITION UPDATES *******************/
Kurt A. O'Hearn
committed
static void Sort_Position_Updates( reax_system * const system, int start, int end,
int dim, mpi_out_data * const out_bufs )
Kurt A. O'Hearn
committed
const reax_atom * const atoms = system->my_atoms;
for ( p = 2 * dim; p < 2 * dim + 2; ++p )
{
out = (rvec*) out_bufs[p].out_atoms;
memcpy( out[i], atoms[ out_bufs[p].index[i] ].x, sizeof(rvec) );
}
/*************** UNPACK POSITION UPDATES ******************/
Kurt A. O'Hearn
committed
static void Unpack_Position_Updates( reax_system * const system, int end, void * const dummy,
int cnt, neighbor_proc * const nbr, int dim )
Kurt A. O'Hearn
committed
int i;
const int start = nbr->atoms_str;
reax_atom * const atoms = system->my_atoms;
Kurt A. O'Hearn
committed
rvec * const src = (rvec*) dummy;
/* adjust coordinates of recved atoms if nbr is a periodic one */
if ( nbr->prdc[dim] )
{
dx = nbr->prdc[dim] * system->big_box.box_norms[dim];
Kurt A. O'Hearn
committed
Kurt A. O'Hearn
committed
int SendRecv( reax_system * const system, mpi_datatypes * const mpi_data, MPI_Datatype type,
int * const nrecv, message_sorter sort_func, unpacker unpack, int clr )
int d, cnt, start, end, max, est_flag, ret;
Kurt A. O'Hearn
committed
mpi_out_data * const out_bufs = mpi_data->out_buffers;
const MPI_Comm comm = mpi_data->comm_mesh3D;
MPI_Request req1, req2;
MPI_Status stat1, stat2;
neighbor_proc *nbr1, *nbr2;
#if defined(DEBUG_FOCUS)
int i, p;
#endif
fprintf( stderr, "p%d sendrecv: entered\n", system->my_rank );
if ( clr == TRUE )
{
Reset_Out_Buffers( mpi_data->out_buffers, system->num_nbrs );
}
est_flag = (mpi_data->in1_buffer == NULL) || (mpi_data->in2_buffer == NULL) ?
TRUE : FALSE;
for ( d = 0; d < 3; ++d )
{
sort_func( system, start, end, d, out_bufs );
start = end;
Kurt A. O'Hearn
committed
nbr1 = &system->my_nbrs[2 * d];
nbr2 = &system->my_nbrs[2 * d + 1];
/* for estimates in1_buffer & in2_buffer will be NULL */
if ( est_flag == TRUE )
Estimate_Init_Storage( system->my_rank, nbr1, nbr2, d,
Kurt A. O'Hearn
committed
&max, nrecv, mpi_data, comm );
#if defined(DEBUG_FOCUS)
for ( p = 2 * d; p < 2 * d + 2; ++p )
{
for ( i = 0; i < out_bufs[p].cnt; ++i )
{
fprintf( stderr, "p%d: out_bufs[%d].index[%d] = %d\n",
system->my_rank, p, i, out_bufs[p].index[i] );
fprintf( stderr, " p%d: atom %6d, x[0] = %10.4f, x[1] = %10.4f, x[2] = %10.4f\n",
system->my_rank,
((boundary_atom *)(out_bufs[p].out_atoms))[i].orig_id,
((boundary_atom *)(out_bufs[p].out_atoms))[i].x[0],
((boundary_atom *)(out_bufs[p].out_atoms))[i].x[1],
((boundary_atom *)(out_bufs[p].out_atoms))[i].x[2] );
}
}
#endif
Kurt A. O'Hearn
committed
ret = MPI_Irecv( mpi_data->in1_buffer, nrecv[2 * d], type, nbr1->rank, 2 * d + 1, comm, &req1 );
Check_MPI_Error( ret, "SendRecv::MPI_Irecv::nbr1" );
Kurt A. O'Hearn
committed
ret = MPI_Irecv( mpi_data->in2_buffer, nrecv[2 * d + 1], type, nbr2->rank, 2 * d, comm, &req2 );
Check_MPI_Error( ret, "SendRecv::MPI_Irecv::nbr2" );
ret = MPI_Send( out_bufs[2 * d].out_atoms, out_bufs[2 * d].cnt, type,
Kurt A. O'Hearn
committed
nbr1->rank, 2 * d, comm );
Check_MPI_Error( ret, "SendRecv::MPI_Send::nbr1" );
ret = MPI_Send( out_bufs[2 * d + 1].out_atoms, out_bufs[2 * d + 1].cnt, type,
Kurt A. O'Hearn
committed
nbr2->rank, 2 * d + 1, comm );
Check_MPI_Error( ret, "SendRecv::MPI_Send::nbr2" );
/* recv and unpack atoms from nbr1 in dimension d */
ret = MPI_Wait( &req1, &stat1 );
Check_MPI_Error( ret, "SendRecv::MPI_Wait::nbr1" );
ret = MPI_Get_count( &stat1, type, &cnt );
Check_MPI_Error( ret, "SendRecv::MPI_Count::nbr1" );
Kurt A. O'Hearn
committed
unpack( system, end, mpi_data->in1_buffer, cnt, nbr1, d );
Kurt A. O'Hearn
committed
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: nbr1: d = %d, end = %d\n", system->my_rank, d, end );
#endif
ret = MPI_Wait( &req2, &stat2 );
Check_MPI_Error( ret, "SendRecv::MPI_Wait::nbr2" );
ret = MPI_Get_count( &stat2, type, &cnt );
Check_MPI_Error( ret, "SendRecv::MPI_Count::nbr2" );
Kurt A. O'Hearn
committed
unpack( system, end, mpi_data->in2_buffer, cnt, nbr2, d );
Kurt A. O'Hearn
committed
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: nbr2: d = %d, end = %d\n", system->my_rank, d, end );
#endif
if ( est_flag == TRUE )
{
system->est_recv = max;
system->est_trans = (max * sizeof(boundary_atom)) / sizeof(mpi_atom);
}
return end;
Kurt A. O'Hearn
committed
void Comm_Atoms( reax_system * const system, control_params * const control,
simulation_data * const data, storage * const workspace,
mpi_datatypes * const mpi_data, int renbr )
real t_start = 0, t_elapsed = 0;
if ( system->my_rank == MASTER_NODE )
if ( renbr == TRUE )
Kurt A. O'Hearn
committed
/* transfer ownership of atoms */
for ( i = 0; i < MAX_NBRS; ++i )
{
nrecv[i] = system->est_trans;
}
system->n = SendRecv( system, mpi_data, mpi_data->mpi_atom_type, nrecv,
Sort_Transfer_Atoms, Unpack_Transfer_Message, TRUE );
Bin_My_Atoms( system, workspace );
Kurt A. O'Hearn
committed
fprintf( stderr, "p%d, step %d: updated local atoms, n=%d\n",
system->my_rank, data->step, system->n );
Kurt A. O'Hearn
committed
/* exchange ghost region info with neighbors */
for ( i = 0; i < MAX_NBRS; ++i )
{
nrecv[i] = system->my_nbrs[i].est_recv;
}
Kurt A. O'Hearn
committed
system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
Sort_Boundary_Atoms, Unpack_Exchange_Message, TRUE );
Kurt A. O'Hearn
committed
fprintf( stderr, "p%d, step %d: exchanged boundary atoms, N=%d\n",
system->my_rank, data->step, system->N );
fprintf( stderr, "p%d: nbr%d(p%d) str=%d cnt=%d end=%d\n",
Kurt A. O'Hearn
committed
system->my_rank, i, system->my_nbrs[i].rank,
system->my_nbrs[i].atoms_str, system->my_nbrs[i].atoms_cnt,
system->my_nbrs[i].atoms_str + system->my_nbrs[i].atoms_cnt );
for ( i = 0; i < MAX_NBRS; ++i )
{
nrecv[i] = system->my_nbrs[i].atoms_cnt;
}
SendRecv( system, mpi_data, mpi_data->mpi_rvec, nrecv,
Sort_Position_Updates, Unpack_Position_Updates, FALSE );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: updated positions\n", system->my_rank );
MPI_Barrier( MPI_COMM_WORLD );
if ( system->my_rank == MASTER_NODE )
{
t_elapsed = Get_Timing_Info( t_start );
data->timing.comm += t_elapsed;
}
Kurt A. O'Hearn
committed
fprintf( stderr, "p%d @ renbr=%d: comm_atoms done\n",
system->my_rank, renbr );
fprintf( stderr, "p%d: system->n = %d, system->N = %d\n",
system->my_rank, system->n, system->N );
//Print_My_Ext_Atoms( system );
//Print_All_GCells( system );
//MPI_Barrier( MPI_COMM_WORLD );