/*
    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/>.
*/
#include "sim.h"

struct S_STATE{
    int n_bonds;
    int n_horizontal_bonds;
    int n_queue;
    int n_traj;
    int end;
};

S_STATE h_state;
int n_grid;
int n_horizontal_grid;

unsigned int *h_grid;
unsigned int *h_horizontal_grid;
unsigned int *h_traj;

#include "random.cu"
#include "kernel.cu"
#include "purge.cu"

void simulate_alloc(dim3 d,bool save_traj)
{
    int s=((d.x*d.y*d.z)*3)>>5;
    s+=2;
    n_grid=s*sizeof(unsigned int);
    h_grid=(unsigned int *)malloc(n_grid);
    memset((void *)h_grid,0,n_grid);
    
    CUDA_SAFE_CALL(cudaMalloc((void **)&s_grid,n_grid));
    CUDA_SAFE_CALL(cudaMemset(s_grid,0,n_grid));
    
    s=(d.x*d.y)>>5;
    s+=2;
    n_horizontal_grid=s*sizeof(unsigned int);
    h_horizontal_grid=(unsigned int *)malloc(n_horizontal_grid);
    memset((void *)h_horizontal_grid,0,n_horizontal_grid);

    CUDA_SAFE_CALL(cudaMalloc((void **)&s_horizontal_grid,n_horizontal_grid));
    CUDA_SAFE_CALL(cudaMemset(s_horizontal_grid,0,n_horizontal_grid));

    CUDA_SAFE_CALL(cudaMalloc((void **)&s_queue,MAX_QUEUE*sizeof(unsigned int)));
    CUDA_SAFE_CALL(cudaMemset(s_queue,0,MAX_QUEUE*sizeof(unsigned int)));
    
    CUDA_SAFE_CALL(cudaMalloc((void **)&s_state,sizeof(S_STATE)));

    if(save_traj){
        CUDA_SAFE_CALL(cudaMalloc((void **)&s_traj,sizeof(unsigned int)*MAX_QUEUE));
        h_traj=(unsigned int *)malloc(sizeof(unsigned int)*MAX_QUEUE);
    }
}

void simulate_save_result(void)
{
    CUDA_SAFE_CALL(cudaMemcpy(h_grid, s_grid, n_grid, cudaMemcpyDeviceToHost));    
    CUDA_SAFE_CALL(cudaMemcpy(h_horizontal_grid, s_horizontal_grid, n_horizontal_grid, cudaMemcpyDeviceToHost));    
}

void simulate_get_state(void)
{
    simulate_get_state_kernel<<<1,1>>>(s_state);
    main_check();
    CUDA_SAFE_CALL(cudaMemcpy(&h_state, s_state, sizeof(S_STATE), cudaMemcpyDeviceToHost));    
}

void simulate_free(bool save_traj)
{
    free((void *)h_grid);
    free((void *)h_horizontal_grid);
    
    CUDA_SAFE_CALL(cudaFree(s_grid));
    CUDA_SAFE_CALL(cudaFree(s_horizontal_grid));
    CUDA_SAFE_CALL(cudaFree(s_queue));
    CUDA_SAFE_CALL(cudaFree(s_state));

    if(save_traj){
        free((void *)h_traj);
        CUDA_SAFE_CALL(cudaFree(s_traj));
    }
}

uint4 simulate_traj_get(int poz)
{
    return make_uint4((h_traj[poz]>>22)&1023,
                      (h_traj[poz]>>12)&1023,
                      (h_traj[poz]>>2)&1023,
                      (h_traj[poz])&3);
}


void simulate_save_traj(FILE *s)
{
    if(!s) return;
    if(h_state.n_traj<1) return;
    CUDA_SAFE_CALL(cudaMemcpy(h_traj, s_traj, h_state.n_traj*sizeof(unsigned int), cudaMemcpyDeviceToHost));
    uint4 a;
    for(int i=0;i<h_state.n_traj;i++){
        a=simulate_traj_get(i);
        fprintf(s,"%d %d %d %d\n",a.x,a.y,a.z,a.w);
    }
}

void simulate_save_results(FILE *s)
{
    if(!s) return;
    fprintf(s,"<RESULTS>\n");
    fprintf(s,"%d %d %d\n",h_state.end,h_state.n_bonds,h_state.n_horizontal_bonds);
}

FILE * simulate_fopen(char *fname,dim3 dimension,float k1,float k2)
{
    FILE *tfile=fopen(fname,"w");
    if(tfile){
        cudaDeviceProp deviceProp;
        cudaGetDeviceProperties(&deviceProp,0);
    
        fprintf(tfile,"<TRHLINA_TRAJECTORY 1.0>\n");
        fprintf(tfile,"%g %g\n",k1,k2);
        fprintf(tfile,"%d %d %d\n",dimension.x,dimension.y,dimension.z);
        fprintf(tfile,"<COMMENT>\nGPU version, Card: %s, RAM: %.1fMB, Warp: %d, Threads per block: %d, Clock rate: %.2fGHz\n",
            deviceProp.name,(float)deviceProp.totalGlobalMem/1024/1024,deviceProp.warpSize,deviceProp.maxThreadsPerBlock,(float)deviceProp.clockRate/1000000);
        fprintf(tfile,"<TRAJECTORY>\n");
    }
    return tfile;
}

void simulate_fclose(FILE *s)
{
    if(!s) return;
    fprintf(s,"<END>\n");
    fclose(s);
}

void simulate_start(char *fname,dim3 dimension,bool save_traj,float k1,float k2,int num_t,int2 num_g)
{
    dim3 gdim;
    FILE *tfile=simulate_fopen(fname,dimension,k1,k2);
    if(!tfile) save_traj=false;
    
    init_gpurand(num_t);
    
    gdim=make_uint3(dimension.x+1,dimension.y+1,dimension.z+1);
    simulate_alloc(gdim,save_traj);
    
    printf("Memory: %.2fMB\n",(double)(n_grid+n_horizontal_grid+MAX_QUEUE*4)/1024/1024);
    
    uint4 now=make_uint4(dimension.x/2,dimension.y/2,dimension.z/2,2);
    simulate_init_kernel<<<1,1>>>(dimension,gdim,save_traj,now,s_queue,s_grid,s_horizontal_grid,s_traj);
    main_check();

    simulate_get_state();
    if(save_traj){
        simulate_save_traj(tfile);
    }
    
    int2 sdim=make_int2(num_t,num_g.x);
    
    int n;
    int purge_period=16;
    for(int count=0;;count++){
        if(h_state.n_queue<1) break;
        n=(h_state.n_queue/(num_t+1));
        
        simulate_kernel<<< make_uint3((n%num_g.x)+1,(n/num_g.x)+1,1), make_uint3(num_t,1,1) >>> 
            (dimension,gdim,sdim,save_traj,k1,k2,s_queue,s_grid,s_horizontal_grid,s_traj,gpurandval);
        main_check();
        
        simulate_get_state();
        
        if(save_traj){
            simulate_save_traj(tfile);
        }
        
        if((h_state.end&15)==15) break;
        
        if((count%purge_period)==0){
            int delta=h_state.n_queue;
            simulate_purge(gdim,num_t);
            delta-=h_state.n_queue;
            if(delta<(num_t<<4)) purge_period+=2;
            if(delta>(num_t<<10)) purge_period-=2;
            if(purge_period<2) purge_period=2;
            if(purge_period>16) purge_period=16;
        }
    }

    simulate_get_state();
    simulate_save_results(tfile);
    printf("Broken bonds: %d\n",h_state.n_bonds);

    simulate_free(save_traj);
    simulate_fclose(tfile);
}
