lattice.cpp

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

	    {
	      //	      send[count]=lat[i][j][constant_index];
	      lat[i][j][constant_index].get_nodes(temp);
	      for(int ii=0; ii<temp.size(); ii++)
		node_send[6*count+ii] = temp[ii];
	      count++;
	    }
      
      //      for(int i=0;i<count;i++)
      //      cout<<count<<" printing send stuff here "<<node_send<<endl;
      
       MPI_Isend(&node_send[0],6*size, 
      		MPI_INT, destination,
      		0,MPI_COMM_WORLD,&req);
           //   MPI_Send(&node_send[0], 6*size, 
//        		MPI_INT, destination,
//        		  0,MPI_COMM_WORLD);
    }
}
//----------------------------------------------------------------
void My_recv(const int &source,const int &size,
	     /* const int &x_len,const int &y_len, const int &z_len,
		const int &constant_index,
		vector< vector< vector <Material> > > &lat,*/
	     MPI_Request &req,// vector<Material> &recv,
	     vector<int> &node_recv)
{

  MPI_Status stat;

  if(source >= 0)
    {
      //      cout <<"source " <<source<<" size "<<6*size<<endl;
      node_recv.resize(6*size);

      MPI_Irecv(&node_recv[0],6*size,
		MPI_INT,source,
		0,MPI_COMM_WORLD,&req);
   //           MPI_Recv(&node_recv[0],6*size,
//        		MPI_INT,source,
//        		0,MPI_COMM_WORLD,&stat);
//              cout<<"here's wha twe got "<< node_recv<<endl;

    }
}
//------------------------------------------------------------------------------------------//
void parse_data(const int &source,const int &size,
	     const int &x_len,const int &y_len, const int &z_len,
	     const int &constant_index,
	     vector< vector< vector <Material> > > &lat,
	     MPI_Request &req, vector<int> &node_recv)
{
  vector<int> temp(6,0);
  int count = 0;
  if(source >= 0)
    {
      //need to save the vector here!
      if(x_len == 0)
  	for(int j=1; j<=y_len;j++)
  	  for(int k=1; k<=z_len; k++)
  	    {
  	      for(int ii=0;ii<temp.size();ii++)
  		temp[ii]=node_recv[ii+6*count];
	      lat[constant_index][j][k].set_nodes(temp);
  	      count++;
  	    }
      else if(y_len == 0)
	for(int i=1; i<= x_len;i++)
	  for(int k=1; k<=z_len; k++)
	    {
	      for(int ii=0;ii<temp.size();ii++)
		temp[ii]=node_recv[ii+6*count];
	      lat[i][constant_index][k].set_nodes(temp);
	      count++;
	    }
      else
	for(int i=1; i<= x_len;i++)
	  for(int j=1; j<=y_len;j++)
	    {
	      for(int ii=0;ii<temp.size();ii++)
		temp[ii]=node_recv[ii+6*count];
	      lat[i][j][constant_index].set_nodes(temp);
	      count++;
	    }
    }
}

//--------------------------------------------------------------------------------------------------------//
double calc_inner_wall( vector< vector< vector< Material > > > &lat,
			const int &startx, const int &endx,
			const int &starty, const int &endy,
			const int &startz, const int &endz,
			const vector <int> &nbrs, int &calc_E_count)
{
  double E_crystal = 0;
  int i,j,k;

  for(i=startx; i<=endx; i++)
    for(j=starty; j<=endy; j++)
      for(k=startz; k<=endz; k++)
	{
	  for(int ii = 0; ii<nbrs.size(); ii++)
	    if(nbrs[ii]==1)
	      {
		E_crystal += lat[i][j][k].calc_E(ii,lat);//calc E for nodes
		calc_E_count+=1;
	      }
	  
	}
  return E_crystal;
}
//-------------------------------------------------------------------------------------------//
void  get_start_and_end(int &startx,int &endx,const vector< vector<int> > &sxy,
			const int &X, const int &x,const int &index)
{
  startx = 1;
  endx = X;

  if(sxy[index][0]==1)//we are at front
    startx = 2;
  if(sxy[index][1]==x)//we are at the end
    endx = X-1;
}
//-------------------------------------------------------------------------------------------//
void  get_start_and_end_interior(int &startx,int &endx,const vector< vector<int> > &sxy,
			const int &X, const int &x,const int &index)
{
  startx = 1;
  endx = X;
  
  if(sxy[index][0]==1)//we are at front
    startx = 2;
  if(sxy[index][1]==x)//we are at the end
    endx = X-1;


  if(index==0)//x-direction
    {
      if(sxy[0][0]==1)
	endx--;//we want to end sooner
      else
	startx++;
    }
//    if(index==1)//y direction
//        {
//          if(sxy[1][0]==1)
//  	  endx--;
//  	else
//  	  startx++;
//        }
  //  if(index==2)//z direction
//      {
//        if(sxy[2][1]!=1)
//  	startx++;
//      }
}
//-------------------------------------------------------------------//
void Lattice::get_W()
{
  int node_total=0, count = 1;
  double K = 8.617E-4;
  int done = 0;
  MPI_Status STAT;
  MPI_Allreduce (&calc_E_count,&node_total, 1, 
  		 MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  w[0][0]=0;//zero out the offsets each time through
  int X = lat.size()-2;
  int Y = lat[0].size()-2;
  int Z = lat[0][0].size()-2;
  for(int i = 1; i <= X; i++)//lattice starts at 1 not 0
    for(int j = 1; j <= Y; j++)
      for(int k = 1; k <= Z; k++)
	{
	  w[count][0] = w[count - 1][0] +
	    exp(-lat[i][j][k].get_delta_E(bottom)/(K*Temperature));//calc E for bottom node
	  w[count][1] = i; w[count][2]=j; w[count][3]=k;w[count][4]=bottom;
   	  count++;
	  
	  w[count][0] = w[count - 1][0] + 
	    exp(-lat[i][j][k].get_delta_E(left  )/(K*Temperature));//calc E for left node
	  w[count][1] = i; w[count][2]=j; w[count][3]=k;w[count][4]=left;
	  count++;
	  
	  w[count][0] = w[count - 1][0] + 
   	    exp(-lat[i][j][k].get_delta_E(front )/(K*Temperature));//calc E for front node
	  w[count][1] = i; w[count][2]=j; w[count][3]=k;w[count][4]=front;
   	  count++;

   	  if(i==x)//then we need to do the right node
   	    {
	      w[count][0] = w[count - 1][0] + exp(-lat[i][j][k].get_delta_E(right)/(K*Temperature));
	      w[count][1] = i; w[count][2]=j; w[count][3]=k;w[count][4]=right;
	      count++;
   	    }
   	  if(j==y)//then we need to do the back
   	    {
  	      w[count][0] = w[count - 1][0] + exp(-lat[i][j][k].get_delta_E(back)/(K*Temperature));
	      w[count][1] = i; w[count][2]=j; w[count][3]=k;w[count][4]=back;
   	      count++;
   	    }
 	  if(k==z)//then we need to do the top node
     	    {
    	      w[count][0] = w[count - 1][0] + exp(-lat[i][j][k].get_delta_E(top)/(K*Temperature));
	      w[count][1] = i; w[count][2]=j; w[count][3]=k;w[count][4]=top;
     	      count++;
     	    }
	}
  if(w.size() != (count - 1))
    {
      w.resize(count-1);
      cout <<"resizing w"<<endl;
    }

  if(myid==0)
    {
      offset[0] = w[count-1][0];
      double temp = offset[0];
       for (int i =1; i<num_procs; i++)//get the offsets from each process
 	{
 	  MPI_Recv(&offset[i],1,MPI_DOUBLE,i,0,MPI_COMM_WORLD,&STAT);
 	  temp += offset[i];
 	  offset[i] = temp;
 	  MPI_Send(&offset[i-1],1,MPI_DOUBLE,i,0,MPI_COMM_WORLD);
 	}
       w_length = w[w.size()-1][0];
       if(num_procs != 1)
	 MPI_Recv(&w_length,1,MPI_DOUBLE,num_procs-1,0,MPI_COMM_WORLD,&STAT);
    }
  else
    { 
      MPI_Send(&w[count-1][0],1,MPI_DOUBLE,0,0,MPI_COMM_WORLD);//send our final value
      MPI_Recv(&w[0][0],1,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&STAT);//store the offset in w[0][0]
      if(myid==(num_procs-1)&(num_procs!=1))
	{
	  double temp=w[0][0]+w[w.size()-1][0];
	  MPI_Send(&temp,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
	}
    }
  MPI_Bcast(&w_length,num_procs, MPI_DOUBLE,0,MPI_COMM_WORLD);


  //lets write the w values to a file so we can see a plot of it
  if(myid==0)
    {
      ofstream w_file ("w.dat",ios::out);//open a file for output
      if (!w_file)
	cout <<"WARNING!!!  couldn't open the energy.dat file"<<endl;
      char temp[1000],*tempstr=temp;
      for (int i = 1; i<w.size();i++)
	w_file << setprecision(15)<< (w[i][0] - w[i-1][0])<<" "<<w[i][1] <<" "<< w[i][2]<<" "<<w[i][3]<<endl;
      w_file.close();
    }

}
//---------------------------------------------------------------------------//
void Lattice::rand_flip(int &index_i, int &index_j, int &index_k, int &proc, int &node)
{
  int middle = (int)(floor((w.size()-1)/2));
  double start = 0;
  double end = w.size()-1;
  int done = 0;
  double pick=0;

  proc =-1;

  if(myid==0)//if we are the root process
    {
      double max_w=offset[offset.size()-1]; //get the maximum w value
      //now we generate a random number to pick which oxygen that we flip through
      srand(time(0));//initialize the random number generator
      pick = rand();//get our random number
      done = 0;
      while(!done)
	{
	  if(pick <= max_w)//this makes the random number <= to max_w
 	    done = 1;
	  else if(max_w < 1E-15) //this is special case
	    {
	      pick = 0;
	      done = 1;
	    }
	  else
	    {
	      int temp = (int)(pick/max_w);
	      pick -= temp*max_w;
	    }
 	}
    }
  MPI_Bcast(&pick,1,MPI_DOUBLE,0,MPI_COMM_WORLD);//broadcast the pick to all the processes
  //  cout <<"random pick = " <<pick <<" MYID =========== " <<myid <<" start = "<<w[0][0]<<" end = "<<(w[0][0]+w[w.size()-1][0])<<endl;
  if(pick==0)
    {
      //just flip any node
      index_i = 1;
      index_j = 1;
      index_k = 1;
      node = 1;
      proc = 0;
      int node_value = lat[index_i][index_j][index_k].get_node(node);
      //now flip the value
      node_value = -node_value;//just flip the sign
      //now we set the node value
      lat[index_i][index_j][index_k].set_node(node_value,node);
    }
  else if((w[0][0] <= pick)&((w[0][0]+ w[w.size()-1][0]) >= pick))//determine who has the node to be flipped
     {
       //       cout <<myid <<" has the oxygen to be flipped  start = "<<w[0][0]<<" end = "<<w[0][0]+w[w.size()-1][0]<<endl;
       //now lets find the actual oxygen to flip
       //use bisection method to get there somewhat fast
       done = 0;
       while(!done)
	{
	  if (pick > (w[0][0]+w[middle][0]))
	    {
	      start = middle;
	      middle = (int)(.5*(start+end));
	      //	      cout <<"too high "<<w[0][0]+w[middle][0]<<" middle = " <<middle <<endl;
	    }
	  else//we are below middle
	    {
	      end = middle;
	      middle = (int)(.5*(start+end));
	      //	      cout <<"too low "<<w[0][0]+w[middle][0]<<" middle = " <<middle<<endl;
	    }
	  if((start == middle)|(end == middle))//then we've found it
	    done = 1;
	}
       middle++;//it picks the one below the one we want
       
       //now we flip the node from its current value
       //      cout <<"middle = "<<middle<<" size = " <<w.size()<<endl;
       index_i =(int)w[middle][1];
       index_j =(int)w[middle][2];
       index_k =(int)w[middle][3];
       node    =(int)w[middle][4];
       //       cout <<"i jk = "<<index_i<<" "<<index_j<<" "<<index_k<<endl;
       cout<<"The w at that node before " <<
   	 exp((-lat[index_i][index_j][index_k].get_delta_E(node))/(Temperature))<<" delta E = " <<
   	 lat[index_i][index_j][index_k].get_delta_E(node) <<endl;
       int node_value = lat[index_i][index_j][index_k].get_node(node);
       //now flip the value
       cout <<"we are going to flip the node from "<<node_value<<endl;
       node_value = -node_value;//just flip the sign
       cout <<"the node is now "<<node_value<<endl;
       //now we set the node value
       lat[index_i][index_j][index_k].set_node(node_value,node);
       proc = myid;
     }

}
//-----------------------------------------------------------------------------------------------------------------------//
void Lattice::adjust_E_crystal(int &index_i, int &index_j, int &index_k, int &proc, int &node)
{
  
  //this is extremely lazy! we can really boost performance by only
  //recomputing the engery of its neighboors and itself
  //its just easier to have it recompute its entire crystal energy

  this->exchange();
  this->inner_comps();
  this->wait();
  this->outer_comps();
  if(myid==proc)
    {
      lat[index_i][index_j][index_k].calc_E(node,lat);
      cout<<"The w at that node after " <<
	exp((-lat[index_i][index_j][index_k].get_delta_E(node))/((8.617E-3)*Temperature))<<" delta E = " <<
	lat[index_i][index_j][index_k].get_delta_E(node) <<endl;
    }
  
}