14 #include "Foundation/console.h"
19 #ifdef __INTEL_COMPILER
30 #include "ntable/src/nt_slab.h"
68 int nx_global,ny_global,nz_global;
69 nx_global=lrint((max[0]-min[0])/range);
70 ny_global=lrint((max[1]-min[1])/range);
71 nz_global=lrint((max[2]-min[2])/range);
72 if ((fabs((
double(nx_global)-(max[0]-min[0])/range)) > 1e-6) ||
73 (fabs((
double(ny_global)-(max[1]-min[1])/range)) > 1e-6) ||
74 (fabs((
double(nz_global)-(max[2]-min[2])/range)) > 1e-6)){
76 std::stringstream msg;
77 msg <<
"size doesn't fit range" << endl;
78 msg <<
"diff x : " << double(nx_global)-(max[0]-min[0])/range << endl;
79 msg <<
"diff y : " << double(ny_global)-(max[1]-min[1])/range << endl;
80 msg <<
"diff z : " << double(nz_global)-(max[2]-min[2])/range << endl;
81 console.
Error() << msg.str() <<
"\n";
82 throw std::runtime_error(msg.str());
85 int nx_min=((nx_global*pcoords[0])/m_comm.
get_dim(0))-1;
86 int ny_min=((ny_global*pcoords[1])/m_comm.
get_dim(1))-1;
87 int nz_min=((nz_global*pcoords[2])/m_comm.
get_dim(2))-1;
89 int nx=(((nx_global*(pcoords[0]+1))/m_comm.
get_dim(0))-nx_min)+1;
90 int ny=(((ny_global*(pcoords[1]+1))/m_comm.
get_dim(1))-ny_min)+1;
91 int nz=(((nz_global*(pcoords[2]+1))/m_comm.
get_dim(2))-nz_min)+1;
93 m_nt=
new NeighborTable<T>(nx,ny,nz,range,alpha,min,nx_min,ny_min,nz_min);
127 int nx_global,ny_global,nz_global;
128 nx_global=lrint((max[0]-min[0])/range);
129 ny_global=lrint((max[1]-min[1])/range);
130 nz_global=lrint((max[2]-min[2])/range);
131 if ((fabs((
double(nx_global)-(max[0]-min[0])/range)) > 1e-6) ||
132 (fabs((
double(ny_global)-(max[1]-min[1])/range)) > 1e-6) ||
133 (fabs((
double(nz_global)-(max[2]-min[2])/range)) > 1e-6)){
135 std::stringstream msg;
136 msg <<
"size doesn't fit range" << endl;
137 msg <<
"diff x : " << double(nx_global)-(max[0]-min[0])/range << endl;
138 msg <<
"diff y : " << double(ny_global)-(max[1]-min[1])/range << endl;
139 msg <<
"diff z : " << double(nz_global)-(max[2]-min[2])/range << endl;
140 throw std::runtime_error(msg.str());
143 int nx_min=((nx_global*pcoords[0])/m_comm.
get_dim(0))-1;
144 int ny_min=((ny_global*pcoords[1])/m_comm.
get_dim(1))-1;
145 int nz_min=((nz_global*pcoords[2])/m_comm.
get_dim(2))-1;
147 int nx=(((nx_global*(pcoords[0]+1))/m_comm.
get_dim(0))-nx_min)+1;
148 int ny=(((ny_global*(pcoords[1]+1))/m_comm.
get_dim(1))-ny_min)+1;
149 int nz=(((nz_global*(pcoords[2]+1))/m_comm.
get_dim(2))-nz_min)+1;
151 m_nt=
new NeighborTable<T>(nx,ny,nz,range,alpha,min,nx_min,ny_min,nz_min);
153 m_xshift=(max-min).X();
154 m_yshift=(max-min).Y();
155 m_zshift=(max-min).Z();
157 m_circ_edge_x_down=(circ[0] && (pcoords[0]==0)) ?
true : false ;
158 m_circ_edge_x_up=(circ[0] && (pcoords[0]==m_comm.
get_dim(0)-1)) ?
true : false ;
170 if(m_nt!=NULL)
delete m_nt;
192 for(
typename vector<T>::const_iterator iter=vp.begin();
207 return m_nt->isInInner(pos);
219 return m_nt->ptr_by_id(
id);
231 return m_nt->getNearestPtr(pos);
243 vector<int> pcoords=m_comm.get_coords();
253 bool circ_edge=
false;
256 if(m_comm.get_dim(0)>1){
258 NTSlab<T> up_boundary_slab=m_nt->yz_slab(m_nt->xsize()-1);
259 up_boundary_slab.erase(up_boundary_slab.begin(),up_boundary_slab.end());
260 NTSlab<T> down_boundary_slab=m_nt->yz_slab(0);
261 down_boundary_slab.erase(down_boundary_slab.begin(),down_boundary_slab.end());
267 up_slab=m_nt->yz_slab(m_nt->xsize()-2);
268 down_slab=m_nt->yz_slab(1);
271 if(m_circ_edge_x_up){
275 for(
typename NTSlab<T>::iterator iter=up_slab.begin();
278 buffer.push_back(*iter);
281 for(
typename vector<T>::iterator iter=buffer.begin();
284 iter->setCircular(
Vec3(-1.0*m_xshift,0.0,0.0));
286 m_comm.shift_cont_packed(buffer,*m_nt,0,1,m_exchg_tag);
288 m_comm.shift_cont_packed(up_slab,*m_nt,0,1,m_exchg_tag);
291 if(m_circ_edge_x_down){
295 for(
typename NTSlab<T>::iterator iter=down_slab.begin();
296 iter!=down_slab.end();
298 buffer.push_back(*iter);
301 for(
typename vector<T>::iterator iter=buffer.begin();
304 iter->setCircular(
Vec3(m_xshift,0.0,0.0));
306 m_comm.shift_cont_packed(buffer,*m_nt,0,-1,m_exchg_tag);
308 m_comm.shift_cont_packed(down_slab,*m_nt,0,-1,m_exchg_tag);
312 if(m_comm.get_dim(1)>1){
314 NTSlab<T> up_boundary_slab=m_nt->xz_slab(m_nt->ysize()-1);
315 up_boundary_slab.erase(up_boundary_slab.begin(),up_boundary_slab.end());
316 NTSlab<T> down_boundary_slab=m_nt->xz_slab(0);
317 down_boundary_slab.erase(down_boundary_slab.begin(),down_boundary_slab.end());
323 up_slab=m_nt->xz_slab(m_nt->ysize()-2);
324 down_slab=m_nt->xz_slab(1);
327 m_comm.shift_cont_packed(up_slab,*m_nt,1,1,m_exchg_tag);
329 m_comm.shift_cont_packed(down_slab,*m_nt,1,-1,m_exchg_tag);
331 cout <<
"circular y shift not implemented" << endl;
335 if(m_comm.get_dim(2)>1){
338 NTSlab<T> up_boundary_slab=m_nt->xy_slab(m_nt->zsize()-1);
339 up_boundary_slab.erase(up_boundary_slab.begin(),up_boundary_slab.end());
340 NTSlab<T> down_boundary_slab=m_nt->xy_slab(0);
341 down_boundary_slab.erase(down_boundary_slab.begin(),down_boundary_slab.end());
344 up_slab=m_nt->xy_slab(m_nt->zsize()-2);
345 down_slab=m_nt->xy_slab(1);
348 m_comm.shift_cont_packed(up_slab,*m_nt,2,1,m_exchg_tag);
350 m_comm.shift_cont_packed(down_slab,*m_nt,2,-1,m_exchg_tag);
352 cout <<
"circular x shift not implemented" << endl;
373 if(m_comm.get_dim(0)>1){
376 exchange_single(rdf,wrtf,m_nt->yz_slab(m_nt->xsize()-2),m_nt->yz_slab(0),0,1);
378 exchange_single(rdf,wrtf,m_nt->yz_slab(1),m_nt->yz_slab(m_nt->xsize()-1),0,-1);
381 if(m_comm.get_dim(1)>1){
384 exchange_single(rdf,wrtf,m_nt->xz_slab(m_nt->ysize()-2),m_nt->xz_slab(0),1,1);
386 exchange_single(rdf,wrtf,m_nt->xz_slab(1),m_nt->xz_slab(m_nt->ysize()-1),1,-1);
389 if(m_comm.get_dim(2)>1){
392 exchange_single(rdf,wrtf,m_nt->xy_slab(m_nt->zsize()-2),m_nt->xy_slab(0),2,1);
394 exchange_single(rdf,wrtf,m_nt->xy_slab(1),m_nt->xy_slab(m_nt->zsize()-1),2,-1);
411 NTSlab<T> send_slab,NTSlab<T> recv_slab,
418 for(
typename NTSlab<T>::iterator iter=send_slab.begin();
419 iter!=send_slab.end();
421 send_data.push_back(((*iter).*rdf)());
425 m_comm.shift_cont_packed(send_data,recv_data,dir,dist,m_exchg_tag);
429 if(recv_data.size()==recv_slab.size()){
431 for(
typename NTSlab<T>::iterator iter=recv_slab.begin();
432 iter!=recv_slab.end();
434 ((*iter).*wrtf)(recv_data[count]);
439 cerr <<
"rank = " << m_comm.rank() <<
":ParallelParticleArray<T>::exchange_single ERROR: size mismatch in received data, recv_data.size()!=recv_slab.size() - (" << recv_data.size() <<
"!=" << recv_slab.size() <<
"). dir = " << dir <<
", dist = " << dist << endl;
455 T* pp=m_nt->ptr_by_id(
id);
475 T* pp=m_nt->ptr_by_id(
id);
490 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
493 if(iter->getTag()==tag){
510 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
513 if(iter->getTag()==tag){
531 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
534 if((iter->getTag() & mask) == (tag & mask)){
554 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
557 if((iter->getTag() & mask) == (tag & mask)){
569 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
582 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
599 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
617 NTBlock<T> InnerBlock=m_nt->inner();
618 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
619 iter!=InnerBlock.end();
641 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
644 cont.push_back(((*iter).*rdf)());
650 const NtBlock &ntBlock
652 : m_ntBlock(ntBlock),
653 m_it(m_ntBlock.begin())
655 m_it = m_ntBlock.begin();
656 m_numRemaining = m_ntBlock.size();
662 return (m_numRemaining > 0);
668 Particle &p = (*m_it);
677 return m_numRemaining;
683 return ParticleIterator(m_nt->inner());
697 NTBlock<T> InnerBlock=m_nt->inner();
698 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
699 iter!=InnerBlock.end();
701 cont.push_back(((*iter).*rdf)());
712 template <
typename P>
715 vector<pair<int,P> > res;
717 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
720 res.push_back(make_pair(iter->getID(),((*iter).*rdf)()));
733 template <
typename P>
736 vector<pair<int,P> > res;
738 NTBlock<T> InnerBlock=m_nt->inner();
739 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
740 iter!=InnerBlock.end();
742 res.push_back(make_pair(iter->getID(),((*iter).*rdf)()));
767 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
770 if((iter->getTag() | mask )==(tag | mask)){
771 cont.push_back(((*iter).*rdf)());
789 NTBlock<T> InnerBlock=m_nt->inner();
790 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
791 iter!=InnerBlock.end();
793 int ptag=iter->getTag();
794 if((ptag & mask )==(tag & mask)){
795 cont.push_back(((*iter).*rdf)());
809 template <
typename P>
812 vector<pair<int,P> > res;
814 for(
typename NeighborTable<T>::iterator iter=m_nt->begin();
817 int id=iter->getID();
818 int ptag=iter->getTag();
819 if(( ptag & mask )==(tag & mask)){
820 res.push_back(make_pair(
id,((*iter).*rdf)()));
837 template <
typename P>
840 vector<pair<int,P> > res;
842 NTBlock<T> InnerBlock=m_nt->inner();
843 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
844 iter!=InnerBlock.end();
846 int id=iter->getID();
847 int ptag=iter->getTag();
848 if((ptag & mask )==(tag & mask)){
849 res.push_back(make_pair(
id,((*iter).*rdf)()));
871 template <
typename P>
873 double dx,
double dy,
double dz,
int nx,
int ny,
int nz)
875 console.
Debug() <<
"forPointsGetNearest" <<
"\n";
880 for(
int ix=0;ix<nx;ix++){
881 for(
int iy=0;iy<ny;iy++){
882 for(
int iz=0;iz<nz;iz++){
883 Vec3 p=orig+double(ix)*dx*u_x+double(iy)*dy*u_y+double(iz)*dz*u_z;
884 cont.push_back(((*(m_nt->getNearestPtr(p))).*rdf)());
889 console.
Debug() <<
"end forPointsGetNearest" <<
"\n";
902 NTSlab<T> slab,slab2;
907 slab=m_nt->yz_slab(m_nt->xsize()-1);
908 slab2=m_nt->yz_slab(m_nt->xsize()-2);
910 slab=m_nt->yz_slab(0);
911 slab2=m_nt->yz_slab(1);
916 slab=m_nt->xz_slab(m_nt->ysize()-1);
917 slab2=m_nt->xz_slab(m_nt->ysize()-2);
919 slab=m_nt->xz_slab(0);
920 slab2=m_nt->xz_slab(1);
925 slab=m_nt->xy_slab(m_nt->zsize()-1);
926 slab2=m_nt->xy_slab(m_nt->zsize()-2);
928 slab=m_nt->xy_slab(0);
929 slab2=m_nt->xy_slab(1);
933 cout <<
"Error: wrong direction " << dir <<
" in getBoundarySlabIds" << endl;
936 for(
typename NTSlab<T>::iterator iter=slab.begin();
939 res.insert(iter->getID());
941 for(
typename NTSlab<T>::iterator iter=slab2.begin();
944 res.insert(iter->getID());
965 slab=m_nt->yz_slab(m_nt->xsize()-2);
967 slab=m_nt->yz_slab(1);
972 slab=m_nt->xz_slab(m_nt->ysize()-2);
974 slab=m_nt->xz_slab(1);
979 slab=m_nt->xy_slab(m_nt->zsize()-2);
981 slab=m_nt->xy_slab(1);
985 cout <<
"Error: wrong direction " << dir <<
" in get2ndSlabIds" << endl;
988 for(
typename NTSlab<T>::iterator iter=slab.begin();
991 res.insert(iter->getID());
1002 template<
typename T>
1005 NTBlock<T> InnerBlock=m_nt->inner();
1006 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
1007 iter!=InnerBlock.end();
1009 pv.push_back(*iter);
1019 template<
typename T>
1022 console.
Debug() <<
"ParallelParticleArray<T>::saveCheckPointData\n";
1025 ost << m_nt->base_point() <<
"\n";
1026 ost << m_nt->base_idx_x() <<
" " << m_nt->base_idx_y() <<
" " << m_nt->base_idx_z() <<
"\n";
1029 ost << getInnerSize() <<
"\n";
1032 NTBlock<T> InnerBlock=m_nt->inner();
1033 for(
typename NTBlock<T>::iterator iter=InnerBlock.begin();
1034 iter!=InnerBlock.end();
1036 iter->saveCheckPointData(ost);
1046 template<
typename T>
1049 console.
Debug() <<
"ParallelParticleArray<T>::loadCheckPointData\n";
1060 if((bp!=m_nt->base_point()) ||
1061 (bix!=m_nt->base_idx_x()) ||
1062 (biy!=m_nt->base_idx_y()) ||
1063 (biz!=m_nt->base_idx_z())){
1064 std::cerr <<
"local area data inconsistet: bp " << bp <<
" / " << m_nt->base_point()
1065 <<
" bix: " << bix <<
" / " << m_nt->base_idx_x()
1066 <<
" biy: " << biy <<
" / " << m_nt->base_idx_y()
1067 <<
" bix: " << biz <<
" / " << m_nt->base_idx_z() << std::endl;
1074 for(
int i=0;i!=nparts;i++){
1076 p.loadCheckPointData(ist);
1077 if(!
isInInner(p.getPos())) std::cerr <<
"not in inner: " << p.getPos() << std::endl;
1089 template<
typename T>
1090 ostream& operator<<(ostream& ost, const ParallelParticleArray<T> &ppa)
1092 ost <<
"--ParallelParticleArray--" << endl;
1093 ost << *(ppa.m_nt) << endl;