lattice.cpp

用monte carlo方法计算化学反应的代码Oxygen Ion Ordering

#include "lattice.h"
#include <stdio.h>
Lattice::Lattice(const int &X,const int &Y,const int &Z, 
		 const int &NUM_PROCS, const int &MYID,
		 vector<int> Dims):
  x(X), y(Y), z(Z),//get the supernode size of the lattice
  num_procs(NUM_PROCS), myid(MYID),//get the number of processors and processor id
  coords(NUM_PROCS), nbrs(NUM_PROCS),//initialize the coordinates and nbr vectors
  dims(3), send(req_size),recv(req_size), 
  node_send(req_size), node_recv(req_size), E_crystal(0),num_procs_bad(num_procs),
  calc_E_count(0),offset(num_procs), w_length(0)
  
{
  decompose(num_procs,dims); //decompose the lattice
  
  for(int i=0;i<num_procs;i++)//allocate the memory for the vectors
    {
      coords[i].resize(3);//x, y, z
      nbrs[i].resize(6);  //top, bottom , left, right, front, back
      //this is looking along the dir of y, with x pointing to the right and z pointing up
    }
  get_cart(dims,coords);//get our cartesian coords
  find_neighboors(nbrs,myid,coords,dims);  //now we need to find the neighbors
  
  //setup memory for our lattice in each supernode

  int start,end;

  sxy.resize(3);//x y z
  sxy[0].resize(2);//x start end 
  sxy[1].resize(2);//y start end 
  sxy[2].resize(2);//z start end 
  MPE_DECOMP1D(x,dims[0],coords[myid][0],sxy[0][0],sxy[0][1]);
  MPE_DECOMP1D(y,dims[1],coords[myid][1],sxy[1][0],sxy[1][1]);
  MPE_DECOMP1D(z,dims[2],coords[myid][2],sxy[2][0],sxy[2][1]);
  cout <<myid<<" coords = " <<coords[myid]<<" start end  ="<<sxy[0]<<"," <<sxy[1]<<" " <<sxy[2]<<endl;
  //check to see if we have a feasible amount of processors
  num_procs_bad=0;
  if(((sxy[0][1] - sxy[0][0]) < 0) | ((sxy[1][1]-sxy[1][0])<0) | ((sxy[2][1]-sxy[2][0])<0))
    {
      cout <<myid<<" BAD NUMER OF PROCESSORS!!!"<<endl;
      num_procs_bad = 1;
    }

  if(myid==0)
    {
      int temp=0;
      for(int i=1; i<num_procs;i++)
	{
	  MPI_Recv(&temp,1,MPI_INT,i,0,MPI_COMM_WORLD,&status[0]);
	  if(temp == 1)
	    num_procs_bad = 1;
	}
    }
  else
    MPI_Send(&num_procs_bad,1,MPI_INT,0,0,MPI_COMM_WORLD);

  MPI_Bcast(&num_procs_bad,1,MPI_INT,0,MPI_COMM_WORLD);

  lat.resize(sxy[0][1]-sxy[0][0] + 3);//resize the x-dir
  for(int i=0;i < lat.size();i++)
    {
      lat[i].resize(sxy[1][1]-sxy[1][0] + 3);//resize the y-dir
      for(int j=0;j < lat[i].size();j++)
	lat[i][j].resize(sxy[2][1]-sxy[2][0] + 3);//resize the z-dir
    }
  
  for(int i=0; i<req_size; i++)
    {
      //send[i].resize((sxy[2][1]-sxy[2][0]+1)*(sxy[2][1]-sxy[2][0]+1));//x*y = y*z = x*z  we should check this
      //      recv[i].resize(send[i].size());//but we will use a cube so its doesn't matter
      //      node_recv[i].resize(6*send[i].size());
    }
  //  cout <<"wall vectors of length " <<send[0].size()<<endl;
  
  
  int X = lat.size()-2;
  int Y = lat[0].size()-2;
  int Z = lat[0][0].size()-2;
  
  int size = 8*X*Y*Z+1;//this is an approximate guess (it should overestimate the size we need this is crappy progammin!!)
  w.resize(size);
  for(int i =0; i<w.size(); i++)
    w[i].resize(5);//this part will hold the index value and node value
  w[0][0] = 0;
};


Lattice::~Lattice()
{;};

//------------------Init lattice -------------------------------------//
int Lattice::Init()
{
  double num_o = 2.5;
  int temp=0;

  lat[0][0][0].Init(lat,num_o);
  return num_procs_bad;
};

//-------------------------------exchange data--------------------
void Lattice::exchange()
{
  int x = sxy[0][1] - sxy[0][0]+1;
  int y = sxy[1][1] - sxy[1][0]+1;
  int z = sxy[2][1] - sxy[2][0]+1;

  for (int i = 0; i<req_size; i++)
    request[i] = MPI_REQUEST_NULL;//set the request to null

  My_send(nbrs[myid][bottom],x*y,x,y,0,1,lat,request[0],node_send[0]);//bottom
  My_send(nbrs[myid][top]   ,x*y,x,y,0,z,lat,request[1],node_send[1]);//top
  My_send(nbrs[myid][left]  ,y*z,0,y,z,1,lat,request[2],node_send[2]);//left
  My_send(nbrs[myid][right] ,y*z,0,y,z,x,lat,request[3],node_send[3]);//right
  My_send(nbrs[myid][front] ,x*z,x,0,z,1,lat,request[4],node_send[4]);//front
  My_send(nbrs[myid][back]  ,x*z,x,0,z,y,lat,request[5],node_send[5]);//back

  My_recv(nbrs[myid][bottom],x*y,request[6],node_recv[0]);//recv from the bottom
  My_recv(nbrs[myid][top]   ,x*y,request[7],node_recv[1]);//recv from the top
  My_recv(nbrs[myid][left]  ,y*z,request[8],node_recv[2]);//recv from the left
  My_recv(nbrs[myid][right] ,y*z,request[9],node_recv[3]);//recv from the right
  My_recv(nbrs[myid][front] ,x*z,request[10],node_recv[4]);//recv from the front
  My_recv(nbrs[myid][back]  ,x*z,request[11],node_recv[5]);//recv from the back

};
//------------------------------------------------------------------------//
void Lattice::wait()
{
  int X = sxy[0][1] - sxy[0][0]+1;
  int Y = sxy[1][1] - sxy[1][0]+1;
  int Z = sxy[2][1] - sxy[2][0]+1;
  MPI_Waitall(req_size,request,status); //wait for our outer edge data
  //now we need to parse the data back into the material class
  parse_data(nbrs[myid][bottom],X*Y,X,Y,0,0  ,lat,request[6],node_recv[0]);//recv from the bottom
  parse_data(nbrs[myid][top]   ,X*Y,X,Y,0,Z+1,lat,request[7],node_recv[1]);//recv from the top
  parse_data(nbrs[myid][left]  ,Y*Z,0,Y,Z,0  ,lat,request[8],node_recv[2]);//recv from the left
  parse_data(nbrs[myid][right] ,Y*Z,0,Y,Z,X+1,lat,request[9],node_recv[3]);//recv from the right
  parse_data(nbrs[myid][front] ,X*Z,X,0,Z,0  ,lat,request[10],node_recv[4]);//recv from the front
  parse_data(nbrs[myid][back]  ,X*Z,X,0,Z,Y+1,lat,request[11],node_recv[5]);//recv from the back
 
};
//-----------------------------------------------------------------------//
void Lattice::inner_comps()
{
  //we go through each node and get the E's 
  //this should be the material class function

  E_crystal = 0;//inner_comps should always be called first so this should be ok

  int X = lat.size()-2;
  int Y = lat[0].size()-2;
  int Z = lat[0][0].size()-2;

  for(int i = 2; i < X; i++)//lattice starts at 1 not 0
    for(int j = 2; j < Y; j++)
      for(int k = 2; k < Z; k++)
	{
	  E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
	  E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
	  E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
	  calc_E_count+=3;
	}
}

//-------------------------------------------------------------------------------//
void Lattice::outer_comps()

{
  //we go through each outer wall node and get the E's
  //again this should be the material class function

  //first do the exterior walls
  this->exterior_walls();
  
  //now lets calc the interior walls (facing eachother)
  this->interior_walls();
  
}
//------------------------------------------------------------------------------------//
void Lattice::interior_walls()
{
  int i=1, j=1, k=1;
  int X = lat.size()-2;
  int Y = lat[0].size()-2;
  int Z = lat[0][0].size()-2;
  int startx= 2, endx =X-1;
  int starty= 2, endy =Y-1;
  int startz= 2, endz =Z-1;
  
  vector<int> nbrs(6,0);
  nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1; //bottom, left and fron
  //the decomposition first splits x then y then z
  //so lets go over the x split first (left and right)
  //left wall center
  if ((sxy[0][0]!=1)&(sxy[0][0]!=x))//as long as we aren't on the exterior
    {
      i=1;   
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,startz,endz,nbrs,calc_E_count);
    }
  if((sxy[0][1]!=x)&(sxy[0][0]!=sxy[0][1]))//right wall, we aren't on exterior and we are longer than 1 column
    {
      i=X;
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,startz,endz,nbrs,calc_E_count);
    }
  
  //now we go through the y direction (front and back walls)
  if ((sxy[1][0]!=1)&(sxy[1][0]!=y))//front wall, as long as we aren't on the exterior
    {
      j=1;   
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,startz,endz,nbrs,calc_E_count);
    }
  if((sxy[1][1]!=y)&(sxy[1][0]!=sxy[1][1]))//back wall, we aren't on exterior and we are longer than 1 column
    {
      j=Y;
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,startz,endz,nbrs,calc_E_count);
    }
  
  //  //no we go through the z direction (top and bottom walls)
  if ((sxy[2][0]!=1)&(sxy[2][0]!=z))//bottm wall, as long as we aren't on the exterior
    {
      k=1;   
      E_crystal += calc_inner_wall(lat,startx,endx,starty,endy,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][1]!=z)&(sxy[2][0]!=sxy[2][1]))//back wall, we aren't on exterior and we are longer than 1 column
    {
      k=Z;
      E_crystal += calc_inner_wall(lat,startx,endx,starty,endy,k,k,nbrs,calc_E_count);
    }
  
  //now do the 4 vertical edges
  if((sxy[0][0] != 1)&(sxy[1][0] != 1))//check the 0,0,0 corner first along x direction
    {
      i = 1;      j = 1;
      E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
    }
  if((sxy[1][0] != 1)&(sxy[0][1] != x))//check the x,0,0 corner 
    {
      i = X;     j  = 1;
      E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
    }  
  if(sxy[0][0]!=sxy[0][1])//we aren't one column
    {
      if((sxy[0][1] != x)&(sxy[1][1] != y))//check the x,y,0 corner
	{
	  i = X;      j = Y;
	  E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
	}
      if((sxy[0][0] != 1)&(sxy[1][1] != y))//check the 0,Y,0 corner 
	{
	  i = 1;      j = Y;
	  E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
	}
    }
  
  
  //now do the upper 4 horizontal edges
  if(sxy[2][1]!=z)
    {
      if(sxy[1][0] != 1)//check the 0,0,0 corner first along x direction
	{
	  j = 1;      k = Z;
	  E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
	}
      if(sxy[0][1] != x)//check the x,0,0 corner 
  	{
  	  i = X;     k  = Z;
  	  E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
  	}  
      if(sxy[1][1] != y)//check the x,y,0 corner
    	{
    	  j = Y;      k = Z;
    	  E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
    	}
      if(sxy[0][0] != 1)//check the 0,Y,0 corner 
    	{
    	  i = 1;      k = Z;
    	  E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
    	}
    }
  //now do the lower 4 horizontal edges
  if(sxy[2][0]!=1)
    {
      if(sxy[1][0] != 1)//check the 0,0,0 corner first along x direction
	{
	  j = 1;      k = 1;
	  E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
	}
      if(sxy[0][1] != x)//check the x,0,0 corner 
  	{
  	  i = X;     k  = 1;
  	  E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
  	}  
      if(sxy[1][1] != y)//check the x,y,0 corner
    	{
    	  j = Y;      k = 1;
    	  E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
    	}
      if(sxy[0][0] != 1)//check the 0,Y,0 corner 
    	{
    	  i = 1;      k = 1;
    	  E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
    	}
    }
  
  //now we do the upper corners 
  if(sxy[2][1]!=z)
    {
      if((sxy[0][0] != 1)&(sxy[1][0] != 1))//check the 0,0,z corner 
	E_crystal += calc_inner_wall(lat,1,1,1,1,Z,Z,nbrs,calc_E_count);
      if((sxy[1][0] != 1)&(sxy[0][1] != x))//check the x,0,z corner 
	E_crystal += calc_inner_wall(lat,X,X,1,1,Z,Z,nbrs,calc_E_count);
      if((sxy[0][1] != x)&(sxy[1][1] != y))//check the x,y,z corner
	E_crystal += calc_inner_wall(lat,X,X,Y,Y,Z,Z,nbrs,calc_E_count);
      if((sxy[0][0] != 1)&(sxy[1][1] != y))//check the 0,Y,z corner 
	E_crystal += calc_inner_wall(lat,1,1,Y,Y,Z,Z,nbrs,calc_E_count);
    }
  //now we do the lower corners 
  if(sxy[2][0]!=1)
    {
      if((sxy[0][0] != 1)&(sxy[1][0] != 1))//check the 0,0,z corner 
	E_crystal += calc_inner_wall(lat,1,1,1,1,1,1,nbrs,calc_E_count);
      if((sxy[1][0] != 1)&(sxy[0][1] != x))//check the x,0,z corner 
	E_crystal += calc_inner_wall(lat,X,X,1,1,1,1,nbrs,calc_E_count);
      if((sxy[0][1] != x)&(sxy[1][1] != y))//check the x,y,z corner
	E_crystal += calc_inner_wall(lat,X,X,Y,Y,1,1,nbrs,calc_E_count);
      if((sxy[0][0] != 1)&(sxy[1][1] != y))//check the 0,Y,z corner 
	E_crystal += calc_inner_wall(lat,1,1,Y,Y,1,1,nbrs,calc_E_count);
    }
  
}
//------------------------------------------------------------------------------------//
void Lattice::exterior_walls()
{
  int i=1, j=1, k=1;
  int X = lat.size()-2;
  int Y = lat[0].size()-2;
  int Z = lat[0][0].size()-2;
  //  vector<int> nbrs(6,0);
  int startx, endx, starty, endy, startz, endz;

  
  //do the front wall sans corners
  if(sxy[1][0]==1)
    {
      vector<int> nbrs(6,0);
      j=1;
      get_start_and_end(startx,endx,sxy,X,x,0);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1; //bottom, left and front
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,startz,endz,nbrs,calc_E_count);
    }

  //back wall sans corners
  if(sxy[1][1]==y)
    {
      vector<int> nbrs(6,0);
      j=Y;
      get_start_and_end(startx,endx,sxy,X,x,0);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1; nbrs[back]=1;//bottom, left,front,back
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,startz,endz,nbrs,calc_E_count);
    }
  
  //left wall sans corners
  if(sxy[0][0]==1)
    {
      vector<int> nbrs(6,0);
      i=1;
      get_start_and_end(starty,endy,sxy,Y,y,1);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;//bottom, left,front