calc_oil_domain.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <cmath>
#include <cassert>

struct set_XY{
  double q[2];
  double normq;
  double sk;
};

class PredByNormq{
public:
  bool operator()(const set_XY& riL,const set_XY& riR) const {
    return riL.normq < riR.normq;
  }
};

class FlicSpc{
private:
  struct double3{
    double x,y,z;
  };
  
  std::vector<double3>  pos;
  double		*rho;
  set_XY		*set_qrho_q;
  set_XY                *set_qrho_qmean;
  int			 w_N,a_N,sys_size;
  int			 Nx,Ny,Nxy;
  int			 all_time_step,time_step;
  double		 sc_L[3],grid_leng;

  std::ofstream fout1,foutmean;
  std::ifstream fin;
  const char*	cur_dir;
  const int     begtime;

  double CalcSk(double qx,double qy){
    double sum_sin=0.0; double sum_cos=0.0;
    for(size_t i=0; i<pos.size(); i++){
      const double the = qx*pos[i].x + qy*pos[i].y;
      sum_cos += cos(the); sum_sin += sin(the);
    }
    const double ret = (sum_cos*sum_cos + sum_sin*sum_sin) / pos.size();
    return ret;
  }

public:
  FlicSpc(const char* cur_dir_,int begtime_):cur_dir(cur_dir_),begtime(begtime_){};
  ~FlicSpc(){
    delete [] rho;
    delete [] set_qrho_q;
    delete [] set_qrho_qmean;
  }
  int ret_all_time_step()const{return all_time_step;}
  int ret_time_step()const{return time_step;}
  void f_manag(){
    std::string ss[5];
    std::string s_cur_dir = cur_dir;
    double buffer_[4];
    
    ss[0] = s_cur_dir + "/macro_data.txt";
    std::ifstream fin1(ss[0].c_str());
    fin1 >> w_N >> a_N >> sc_L[0] >> buffer_[0] >> all_time_step >> time_step;
    sc_L[1] = sc_L[0];
    sc_L[2] = sc_L[0];
    sys_size = w_N + a_N;
    
    ss[1] = s_cur_dir + "/macro_IF.txt";
    std::ifstream fin2(ss[1].c_str());
    double tempera = 0.0;
    fin2 >> buffer_[0] >> tempera >> buffer_[1] >>  buffer_[2] >> buffer_[3] >> grid_leng;
    
    ss[2] = s_cur_dir + "/particle_data.txt";
    fin.open(ss[2].c_str());
    
    ss[3] = s_cur_dir + "/spectrum.txt";
    fout1.open(ss[3].c_str());
    
    ss[4] = s_cur_dir + "/meanspect.txt";
    foutmean.open(ss[4].c_str());
  }

  void SetParam(){
    grid_leng *= 2.0;
    Nx = static_cast<int>(sc_L[0]/grid_leng);
    Ny = Nx;
    Nxy = Nx*Ny;
    grid_leng = sc_L[0]/Nx;
  }

  void  LoadPosition(){
    double buf_pos[3]; int pprop; int chem;
    for(int i=0; i<sys_size; i++){
      //NOTE: this is valid only when membrane plane is parallel to XZ plane.
      fin >> buf_pos[0] >> buf_pos[2] >> buf_pos[1] >> pprop >> chem;
    
      //NOTE: this is valid only when membrane plane is parallel to XY plane.
      //fin >> buf_pos[0] >> buf_pos[1] >> buf_pos[2] >> pprop >> chem;
      
      if(pprop == 2 && chem == false){
	const double3 temp_pos = {buf_pos[0],buf_pos[1],buf_pos[2]};
	pos.push_back(temp_pos);
      }
    }
  }

  void SortbyNormq(int curtime){
    for(int ix=0; ix<Nx; ix++){
      for(int iy=0; iy<Ny/2+1; iy++){
	//      const int		hash	   = ix + Nx*iy;
	const int         hash       = iy + (Ny/2+1)*ix;
	const double	qx	   = (ix < Nx/2 + 1) ? 2.*M_PI*ix/(Nx*grid_leng) : 2.*M_PI*(ix-Nx)/(Nx*grid_leng);
	const double	qy	   = 2.*M_PI*iy/(Nx*grid_leng);
	const double	norm_q	   = sqrt(qx*qx + qy*qy);
	set_qrho_q[hash].q[0]        = qx;
	set_qrho_q[hash].q[1]        = qy;
      
	set_qrho_q[hash].normq	   = norm_q;
	set_qrho_q[hash].sk	   = CalcSk(qx,qy);
      
	if(curtime > begtime) set_qrho_qmean[hash].sk   += set_qrho_q[hash].sk;
      }
    }
    std::sort(set_qrho_q,set_qrho_q+(Ny/2+1)*Nx,PredByNormq());
  }
  void MemAlloc(void){
    rho            = new double [Nxy];
    set_qrho_q	   = new set_XY [(Ny/2+1)*Nx];
    set_qrho_qmean = new set_XY [(Ny/2+1)*Nx];

    for(int ix=0; ix<Nx; ix++){
      for(int iy=0; iy<Ny/2+1; iy++){
	const int       hash       = iy + (Ny/2+1)*ix;
	const double	qx	   = (ix < Nx/2 + 1) ? 2. * M_PI*ix / (Nx*grid_leng) : 2.*M_PI*(ix-Nx)/(Nx*grid_leng);
	const double	qy	   = 2.*M_PI*iy/(Nx*grid_leng);
	const double	norm_q	   = sqrt(qx*qx + qy*qy);
	set_qrho_q[hash].q[0]        = qx;
	set_qrho_q[hash].q[1]        = qy;
	set_qrho_q[hash].normq = norm_q;

	set_qrho_qmean[hash].q[0]     = qx;
	set_qrho_qmean[hash].q[1]     = qy;
	set_qrho_qmean[hash].normq    = norm_q;
	set_qrho_qmean[hash].sk   = 0.;
      }
    }
    std::sort(set_qrho_q,set_qrho_q + (Ny/2 + 1) * Nx,PredByNormq());
  }
  void PrintQnorm(void){
    for(int i=0; i<Nx*(Ny/2 + 1); i++){
      fout1 << set_qrho_q[i].normq << " ";
    }
    fout1 << std::endl;
  }
  void AllClear(void){
    for(int i=0; i<Nxy; i++){
      rho[i] = 0.0;
    }
    pos.clear();
  }
  int GenHash(const double* a){
    int indx = static_cast<int>(a[0]/grid_leng);
    int indy = static_cast<int>(a[1]/grid_leng);
    if(indx == Nx) indx--;
    if(indy == Ny) indy--;
    const int hash = Ny*indx + indy;
    assert(hash>=0 && hash<Nxy);
    return hash;
  }
  void WrtFileDat(void){
    for(int i=0; i< Nx*(Ny/2+1); i++){
      fout1 << set_qrho_q[i].sk << " ";
    }
    fout1 << std::endl;
  }
  void TakeMeanTime(void){
    int divis = (begtime/time_step)*time_step;
    divis += time_step;
    int count = (all_time_step - divis)/time_step;
    for(int i=0; i<Nx*(Ny/2+1); i++){
      set_qrho_qmean[i].sk /= count;
    }
  }
  void WrtFileDat_mean(void){
    for(int i=0; i<Nx*(Ny/2+1); i++){
      foutmean << set_qrho_qmean[i].normq << " " << set_qrho_qmean[i].sk << std::endl;
    }
    std::cout << "spectrum element number " << Nx * (Ny/2 + 1) << std::endl;
  }
};

int main(int argc, char *argv[])
{
  if(argc <= 2){
    std::cout << "usage:" << std::endl;
    std::cout << "argv[1] = target directory" << std::endl;
    std::cout << "argv[2] = beg time point" << std::endl;
    std::cout << "argc > 2" << std::endl;
    exit(1);
  }
  
  const int beg_time = atoi(argv[2]);
  const char* cur_dir = argv[1];
  
  FlicSpc Flispc(cur_dir,beg_time);
  
  Flispc.f_manag();
  Flispc.SetParam();
  Flispc.MemAlloc();
  Flispc.PrintQnorm();
  Flispc.AllClear();

  const int	all_time_step = Flispc.ret_all_time_step();
  const int	time_step     = Flispc.ret_time_step();

  for(int time=0; time<all_time_step; time=time+time_step){
    Flispc.AllClear();
    Flispc.LoadPosition();
    Flispc.SortbyNormq(time);
    if(time > beg_time) Flispc.WrtFileDat(); 
  }
   
  Flispc.TakeMeanTime();
  Flispc.WrtFileDat_mean();
   
  std::cout << "./output_dir/spectrum.txt  -> spectrum time evolution"   << std::endl;
  std::cout << "./output_dir/meanspect.txt -> mean value of spectrum "   << std::endl;
  std::cout << "beg time end time " <<  beg_time <<  " " <<  all_time_step  <<  std::endl;
  return 0;
}