/*
    Simulacia sirenia trhliny
    Copyright (C) 2008 Samuel Kupka

    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 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
__device__ S_STATE *s_state;

__device__ int n_bonds,n_horizontal_bonds;
__device__ int n_queue;
__device__ int n_traj;
__device__ int n_state;

__device__ unsigned int *s_traj; // trajectory
__device__ unsigned int *s_queue;
__device__ unsigned int *s_grid;
__device__ unsigned int *s_horizontal_grid;

__constant__ __device__ int s_tab[3][12][4]={
    {
        { 0, 0, 0,1},
        { 0, 0, 0,2},
        { 0,-1, 0,1},
        { 0, 0,-1,2},
        
        {0,0,-1,0},
        {0,0,1,0},
        {0,-1,0,0},
        {0,1,0,0},
        
        {1,0,0,1},
        {1,0,0,2},
        {1,-1,0,1},
        {1,0,-1,2}
    },
    {
        {0,0,0,0},
        {0,0,0,2},
        {-1,0,0,0},
        {0,0,-1,2},
        
        {-1,0,0,1},
        {1,0,0,1},
        {0,0,-1,1},
        {0,0,1,1},
        
        {0,1,0,0},
        {0,1,0,2},
        {-1,1,0,0},
        {0,1,-1,2}
    },
    {
        {0,0,0,0},
        {0,0,0,1},
        {-1,0,0,0},
        {0,-1,0,1},
        
        {-1,0,0,2},
        {1,0,0,2},
        {0,-1,0,2},
        {0,1,0,2},
        
        {0,0,1,0},
        {0,0,1,1},
        {-1,0,1,0},
        {0,-1,1,1}
    }
};

__device__ inline int simulate_check_border(dim3 d,uint4 now)
{
    if(now.x<1||now.y<1||now.z<1||now.x>=d.x||now.y>=d.y||now.z>=d.z){
        return 1;
    }
    return 0;
}

__device__ uint4 simulate_queue_get(int poz,unsigned int *s_queue)
{
    if(poz>=MAX_QUEUE) poz=0;
    uint4 v=make_uint4((s_queue[poz]>>22)&1023,
                      (s_queue[poz]>>12)&1023,
                      (s_queue[poz]>>2)&1023,
                      (s_queue[poz])&3);
    return v;                      
}

__device__ void simulate_queue_set(uint4 v,unsigned int *s_queue)
{
    int poz=atomicAdd(&n_queue,1);
    if(poz>=MAX_QUEUE) return;
    s_queue[poz]=(v.x<<22)
                +(v.y<<12)
                +(v.z<<2)
                +v.w;
}

__device__ void simulate_traj_set(uint4 v,unsigned int *s_traj)
{
    int poz=atomicAdd(&n_traj,1);
    s_traj[poz]=(v.x<<22)
                +(v.y<<12)
                +(v.z<<2)
                +v.w;
}

__device__ uint4 simulate_translate(uint4 d,char o)
{
    return make_uint4(d.x+s_tab[d.w][o][0],
                      d.y+s_tab[d.w][o][1],
                      d.z+s_tab[d.w][o][2],
                      s_tab[d.w][o][3]);
}

__device__ uint2 simulate_grid_translate(dim3 d,uint4 p)
{
    unsigned int s=(p.x+(p.y*d.x)
                  +(p.z*d.x*d.y))*3
                  +p.w;
    uint2 v;
    v.x=s>>5;
    v.y=s&31;
    return v;
}

__device__ uint2 simulate_horizontal_grid_translate(dim3 d,uint4 p)
{
    unsigned int s=p.x
                  +p.y*d.x;
    uint2 v;
    v.x=s>>5;
    v.y=s&31;
    return v;
}

__device__ void simulate_set_grid_bit(dim3 d,uint4 p,bool save_traj,unsigned int *s_grid,unsigned int *s_traj)
{
    uint2 v=simulate_grid_translate(d,p);
    unsigned int bit=1<<(v.y);
    if(!(atomicOr(s_grid+v.x,bit)&bit)){
        atomicAdd(&n_bonds,1);
        if(p.x==0) n_state|=1;
        if(p.y==0) n_state|=2;
        if(p.x==d.x-1) n_state|=4;
        if(p.y==d.y-1) n_state|=8;
        if(save_traj){
            simulate_traj_set(p,s_traj);
        }
    }
}

__device__ void simulate_set_horizontal_grid_bit(dim3 d,uint4 p,unsigned int *s_horizontal_grid)
{
    uint2 v=simulate_horizontal_grid_translate(d,p);
    unsigned int bit=1<<(v.y);
    if(!(atomicOr(s_horizontal_grid+v.x,bit)&bit)){
        atomicAdd(&n_horizontal_bonds,1);
    }
}

__device__ int simulate_get_horizontal_grid_bit(dim3 d,uint4 p,unsigned int *s_horizontal_grid)
{
    uint2 v=simulate_horizontal_grid_translate(d,p);
    if(s_horizontal_grid[v.x]&(1<<(v.y))) return 1;
    return 0;
}

__device__ int simulate_get_grid_bit(dim3 d,uint4 p,unsigned int *s_grid)
{
    uint2 v=simulate_grid_translate(d,p);
    return (s_grid[v.x]>>v.y)&1;
}

__global__ void simulate_kernel(dim3 dimension,dim3 gdim,int2 sdim,bool save_traj,float k1,float k2,unsigned int *s_queue,unsigned int *s_grid,unsigned int *s_horizontal_grid,unsigned int *s_traj,unsigned int *gpurandval)
{
    int poz=threadIdx.x+((blockIdx.y*sdim.y)+blockIdx.x)*sdim.x;
    int bit=1;
    int nb,nhb;
    int s;
    uint4 now,next;
    float now_energy,next_energy;

    if(poz<n_queue){
        now=simulate_queue_get(poz,s_queue);
        next=simulate_translate(now,gpurand_kernel(gpurandval)%12);
        bit=simulate_get_grid_bit(gdim,next,s_grid);
        nb=n_bonds;
        nhb=n_horizontal_bonds;
    }
    __syncthreads();
    if(!bit){
        s=0;
        now_energy=(float)(nb)*k1-k2*powf((float)nhb,1.5f);
        if(next.w==2&&simulate_get_horizontal_grid_bit(gdim,next,s_horizontal_grid)==0){
            s=1;
        }
        next_energy=(float)(nb+1)*k1-k2*powf((float)(nhb+s),1.5f);

        char crack=0;
        if(next_energy<now_energy){
            crack=1;
        }
        else{
            if( round(((float)MAX_GPURAND/expf(next_energy-now_energy))) > gpurand_kernel(gpurandval) ){
                crack=1;
            }
        }
        
        if(crack){
            if(s){
                simulate_set_horizontal_grid_bit(gdim,next,s_horizontal_grid);
            }
            
            if(!simulate_check_border(dimension,next)){
                simulate_queue_set(next,s_queue);
            }
            
            simulate_set_grid_bit(gdim,next,save_traj,s_grid,s_traj);
        }
    }
}

__global__ void simulate_init_kernel(dim3 dimension,dim3 gdim,bool save_traj,uint4 now,unsigned int *s_queue,unsigned int *s_grid,unsigned int *s_horizontal_grid,unsigned int *s_traj)
{
    n_queue=0;
    n_bonds=0;
    n_horizontal_bonds=0;
    n_traj=0;
    n_state=0;
    
    if(now.w==2){
        simulate_set_horizontal_grid_bit(gdim,now,s_horizontal_grid);
    }
    simulate_set_grid_bit(gdim,now,save_traj,s_grid,s_traj);
    simulate_queue_set(now,s_queue);
}

__global__ void simulate_get_state_kernel(S_STATE *x)
{
	x[0].n_bonds=n_bonds;
	x[0].n_horizontal_bonds=n_horizontal_bonds;
	x[0].n_queue=n_queue;
	x[0].n_traj=n_traj;
	x[0].end=n_state;
	n_traj=0;
}
