<commit_before>/*!
 * \file    lbmFlowAroundCylinder.cu
 * \brief   GPU (Cuda) and CPU version running the same code for floating point computation debugging...
 * \author  <NAME>
 * \date    22.01.2017
 */

#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#include <libgen.h>
#include <pgm.h>
#include <timing.h>

#define RE       220.0       // Reynolds number
#define NX       420         // Numer of lattice nodes (width)
#define NY       180         // Numer of lattice nodes (height)
#define LY       ((NY) - 1)  // Height of the domain in lattice units
#define CX       ((NX) / 4)  // X coordinates of the cylinder
#define CY       ((NY) / 2)  // Y coordinates of the cylinder
#define R        ((NY) / 9)  // Cylinder radius
#define ULB      0.04        // Velocity in lattice units
#define NULB     ((ULB) * (R) / (RE))   // Viscoscity in lattice units
#define OMEGA    ((double)1. / (3*(NULB)+0.5))  // Relaxation parameter

#define NB_BLOCKS    1
//#define NB_THREADS 100

#define SQUARE(a) ((a)*(a))
#define GPU_SQUARE(a) (__dmul_rn(a,a))
#define INDEX_2D_FROM_1D(x, y, i) do { (y) = (i)/(NX), (x) = (i)%(NX); } while (0)

typedef enum { OUT_NONE, OUT_FIN, OUT_IMG } out_mode;

typedef struct {
    bool obstacles[NX][NY];  // Should reside in lbm_consts but is too big for constant memory
    double u[NX][NY][2];
    double feq[NX][NY][9];
    double fin[NX][NY][9];
    double fout[NX][NY][9];
    double rho[NX][NY];
    double vel[NX][NY][2];
} lbm_vars;

typedef struct {
    size_t col[3][3];
    size_t opp[9];
    ssize_t v[2][9];
    double t[9];
} lbm_consts;

#ifdef COMPUTE_ON_CPU
// Tweak the code to run on CPU
#define cudaMalloc(dst_ptr, size) do { *(dst_ptr) = (lbm_vars*)malloc(size); } while(0)
#define cudaMemcpy(dst, src, size, mode) memcpy(dst, src, size)
#define cudaMemcpyToSymbol(dst, src, size) memcpy(&dst, src, size)
#define cudaFree(ptr) free(ptr)

#define HANDLE_ERROR(ans) ans
#define HANDLE_KERNEL_ERROR(...) do { __VA_ARGS__; } while(0)

#define fory(...) for (int y = 0; y < NY; ++y) { __VA_ARGS__; }
#define forxy(...) fory(for (int x = 0; x < NX; ++x) { __VA_ARGS__; })

#define RUN_KERNEL(kernel, th1, th2, ...) forxy(kernel(__VA_ARGS__, x, y))

#else
// Code for GPU usage only
#define HANDLE_ERROR(ans) (handleError((ans), __FILE__, __LINE__))
inline void handleError(cudaError_t code, const char *file, int line)
{
   if (code != cudaSuccess) {
      fprintf(stderr,"CUDA assert: %s %s %d\n", cudaGetErrorString(code), file, line);
      exit(EXIT_FAILURE);
   }
}

#define HANDLE_KERNEL_ERROR(...) \
do {                                         \
    __VA_ARGS__;                             \
    HANDLE_ERROR( cudaPeekAtLastError() );   \
    HANDLE_ERROR( cudaDeviceSynchronize() ); \
} while(0)

#define RUN_KERNEL(kernel, th1, th2, ...) HANDLE_KERNEL_ERROR( kernel<<<dimBlock, dimGrid>>>(__VA_ARGS__) )
#endif

// Constants

#ifndef COMPUTE_ON_CPU
__constant__ 
#endif
lbm_consts d_consts;

const ssize_t V[][9] = {
    { 1, 1, 1, 0, 0, 0,-1,-1,-1 },
    { 1, 0,-1, 1, 0,-1, 1, 0,-1 }
};
const double T[] = { 1./36, 1./9, 1./36, 1./9, 4./9, 1./9, 1./36, 1./9, 1./36 };


/**
 * Setup: cylindrical obstacle and velocity inlet with perturbation
 * Creation of a mask with boolean values, defining the shape of the obstacle.
 */
static void initObstacles(lbm_vars* vars)
{
    for (int x = 0; x < NX; x++) {
        for (int y = 0; y < NY; y++) {
            vars->obstacles[x][y] = false;
        }
    }
}

/**
 * Initial velocity profile: almost zero, with a slight perturbation to trigger
 * the instability.
 */
static void initVelocity(lbm_vars* vars)
{
    for (int d = 0; d < 2; d++) {
        for (int x = 0; x < NX; x++) {
            for (int y = 0; y < NY; y++) {
                vars->vel[x][y][d] = (1-d) * ULB * (1 + 0.0001 * sin( y / (double)LY * 2 * M_PI) );
            }
        }
    }
}

static void initRho(lbm_vars* vars)
{
    for (int x = 0; x < NX; x++) {
        for (int y = 0; y < NY; y++) {
            vars->rho[x][y] = 1.0;
        }
    }
    vars->rho[NX/2][NY/2] = 2.0;   
}

static void initCol(size_t* col, ssize_t v0)
{
    for (int f = 0, i = 0; f < 9 && i < 3; f++) {
        if (V[0][f] == v0) {
            col[i++] = f;
        }
    }
}

static void initOpp(size_t* opp)
{
    for (int f = 0; f < 9; f++) {
        for (int g = 0; g < 9; g++) {
            if (V[0][f] == -V[0][g] && V[1][f] == -V[1][g]) {
                opp[f] = g;
                break;
            }
        }
    }
}

#define EQUILIBRIUM_BODY(v, t) \
    do {                                                                 \
        double usqr = 3./2 * ( SQUARE(u[0]) + SQUARE(u[1]) );            \
                                                                         \
        for (int f = 0; f < 9; f++) {                                    \
            double cu = 3 * ( v[0][f] * u[0] + v[1][f] * u[1] );         \
            feq[f] = rho * t[f] * ( 1 + cu + 0.5 * SQUARE(cu) - usqr );  \
        }                                                                \
    } while(0)

#define GPU_EQUILIBRIUM_BODY(v, t) \
    do {                                                                                \
        double usqr = __dmul_rn(3./2, __dadd_rn( GPU_SQUARE(u[0]), GPU_SQUARE(u[1]) )); \
                                                                                        \
        for (int f = 0; f < 9; f++) {                                                   \
            double cu = 3 * ( v[0][f] * u[0] + v[1][f] * u[1] );                        \
            feq[f] = rho * t[f] * ( 1 + cu + 0.5 * SQUARE(cu) - usqr );                 \
        }                                                                               \
    } while(0)


#ifndef COMPUTE_ON_CPU
__host__ 
#endif
static void h_equilibrium(double* feq, double rho, double* u)
{
    EQUILIBRIUM_BODY(V, T);
}

#ifndef COMPUTE_ON_CPU
__device__ 
#endif
static void d_equilibrium(double* feq, double rho, double* u)
{
#ifdef COMPUTE_ON_CPU
    EQUILIBRIUM_BODY(d_consts.v, d_consts.t);
#else
    GPU_EQUILIBRIUM_BODY(d_consts.v, d_consts.t);
#endif
}

#ifndef COMPUTE_ON_CPU
__device__ 
#endif
static void macroscopic(double* fin, double* rho, double* u)
{
    
    *rho = u[0] = u[1] = 0;

    for (int f = 0; f < 9; f++) {
        *rho += fin[f];

        u[0] += d_consts.v[0][f] * fin[f];
        u[1] += d_consts.v[1][f] * fin[f];
    }
    
    u[0] /= *rho;
    u[1] /= *rho;
}

#ifdef COMPUTE_ON_CPU
void lbm_computation(lbm_vars *d_vars, int x, int y) 
{
#else
__global__ void lbm_computation(lbm_vars *d_vars)
{
    for (int idx = threadIdx.x + blockIdx.x * blockDim.x; idx < NX*NY; idx += blockDim.x * gridDim.x) {
        int x, y;
        INDEX_2D_FROM_1D(x, y, idx);
#endif

    // Right wall: outflow condition.
    if (x == NX-1) {
        for (int i = 0; i < 3; i++) {
            int f = d_consts.col[2][i];
            d_vars->fin[NX-1][y][f] = d_vars->fin[NX-2][y][f];
        }
    }

    // Compute macroscopic variables, density and velocity
    macroscopic(d_vars->fin[x][y], &d_vars->rho[x][y], d_vars->u[x][y]);
    
    if (x == 0) {
        // Left wall: inflow condition
        for (size_t d = 0; d < 2; d++) {
            d_vars->u[0][y][d] = d_vars->vel[0][y][d];
        }        

        // Calculate the density
        double s2 = 0, s3 = 0;
        for (size_t i = 0; i < 3; i++) {
            s2 += d_vars->fin[0][y][d_consts.col[1][i]];
            s3 += d_vars->fin[0][y][d_consts.col[2][i]];
        }
        d_vars->rho[0][y] = 1./(1 - d_vars->u[0][y][0]) * (s2 + 2*s3);
    }
   
    // Compute equilibrium
    d_equilibrium(d_vars->feq[x][y], d_vars->rho[x][y], d_vars->u[x][y]);

    if (x == 0) {
        for (size_t i = 0; i < 3; i++) {
            size_t f = d_consts.col[0][i];
            d_vars->fin[0][y][f] = d_vars->feq[0][y][f] + d_vars->fin[0][y][d_consts.opp[f]] - d_vars->feq[0][y][d_consts.opp[f]];
        }
    }
    
    for (size_t f = 0; f < 9; f++) {
        if (d_vars->obstacles[x][y]) {
            // Bounce-back condition for obstacle
            d_vars->fout[x][y][f] = d_vars->fin[x][y][d_consts.opp[f]];
        } else {
            // Collision step
#ifdef COMPUTE_ON_CPU
            d_vars->fout[x][y][f] = d_vars->fin[x][y][f] - OMEGA * (d_vars->fin[x][y][f] - d_vars->feq[x][y][f]);
#else
            d_vars->fout[x][y][f] = __dadd_rn(__dmul_rn(-OMEGA, __dadd_rn(d_vars->fin[x][y][f], - d_vars->feq[x][y][f])), d_vars->fin[x][y][f]);
#endif
        }
    }

#ifndef COMPUTE_ON_CPU
    }
#endif
}

#ifdef COMPUTE_ON_CPU
void lbm_streaming(lbm_vars *d_vars, int x, int y)
{
#else
__global__ void lbm_streaming(lbm_vars *d_vars)
{
    for (int idx = threadIdx.x + blockIdx.x * blockDim.x; idx < NX*NY; idx += blockDim.x * gridDim.x) {
        int x, y;
        INDEX_2D_FROM_1D(x, y, idx);
#endif

    // Streaming step
    for (size_t f = 0; f < 9; f++) {
        ssize_t x_dst = x + d_consts.v[0][f];
        ssize_t y_dst = y + d_consts.v[1][f];
        if (0 <= x_dst && x_dst < NX && 0 <= y_dst && y_dst < NY) {
            d_vars->fin[x_dst][y_dst][f] = d_vars->fout[x][y][f];
        }
    }

#ifndef COMPUTE_ON_CPU
    }
#endif
}

void output_variables(char* filename, double var[NX][NY][9]) 
{
    FILE* file = fopen(filename, "w");
    for (size_t x = 0; x < NX; x++) {
        for (size_t y = 0; y < NY; y++) {
            for (size_t f = 0; f < 9; ++f) {
                fprintf(file, "%64.60f\n", var[x][y][f]);
            }
        }
    }
    fclose(file);
}

void output_image(char* filename, double u[NX][NY][2]) 
{
    pgm_image* pgm = pgm_create(NX, NY);

    for (size_t x = 0; x < NX; x++) {
        for (size_t y = 0; y < NY; y++) {
            double vel = sqrt( SQUARE(u[x][y][0]) + SQUARE(u[x][y][1]) );
            int color =  255 * min(vel * 10, 1.0);<commit_msg>Set velocities to zero and remove right/left walls<commit_after>/*!
 * \file    lbmFlowAroundCylinder.cu
 * \brief   GPU (Cuda) and CPU version running the same code for floating point computation debugging...
 * \author  <NAME>
 * \date    22.01.2017
 */

#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#include <libgen.h>
#include <pgm.h>
#include <timing.h>

#define RE       220.0       // Reynolds number
#define NX       420         // Numer of lattice nodes (width)
#define NY       180         // Numer of lattice nodes (height)
#define LY       ((NY) - 1)  // Height of the domain in lattice units
#define CX       ((NX) / 4)  // X coordinates of the cylinder
#define CY       ((NY) / 2)  // Y coordinates of the cylinder
#define R        ((NY) / 9)  // Cylinder radius
#define ULB      0.04        // Velocity in lattice units
#define NULB     ((ULB) * (R) / (RE))   // Viscoscity in lattice units
#define OMEGA    ((double)1. / (3*(NULB)+0.5))  // Relaxation parameter

#define NB_BLOCKS    1
//#define NB_THREADS 100

#define SQUARE(a) ((a)*(a))
#define GPU_SQUARE(a) (__dmul_rn(a,a))
#define INDEX_2D_FROM_1D(x, y, i) do { (y) = (i)/(NX), (x) = (i)%(NX); } while (0)

typedef enum { OUT_NONE, OUT_FIN, OUT_IMG } out_mode;

typedef struct {
    bool obstacles[NX][NY];  // Should reside in lbm_consts but is too big for constant memory
    double u[NX][NY][2];
    double feq[NX][NY][9];
    double fin[NX][NY][9];
    double fout[NX][NY][9];
    double rho[NX][NY];
    double vel[NX][NY][2];
} lbm_vars;

typedef struct {
    size_t col[3][3];
    size_t opp[9];
    ssize_t v[2][9];
    double t[9];
} lbm_consts;

#ifdef COMPUTE_ON_CPU
// Tweak the code to run on CPU
#define cudaMalloc(dst_ptr, size) do { *(dst_ptr) = (lbm_vars*)malloc(size); } while(0)
#define cudaMemcpy(dst, src, size, mode) memcpy(dst, src, size)
#define cudaMemcpyToSymbol(dst, src, size) memcpy(&dst, src, size)
#define cudaFree(ptr) free(ptr)

#define HANDLE_ERROR(ans) ans
#define HANDLE_KERNEL_ERROR(...) do { __VA_ARGS__; } while(0)

#define fory(...) for (int y = 0; y < NY; ++y) { __VA_ARGS__; }
#define forxy(...) fory(for (int x = 0; x < NX; ++x) { __VA_ARGS__; })

#define RUN_KERNEL(kernel, th1, th2, ...) forxy(kernel(__VA_ARGS__, x, y))

#else
// Code for GPU usage only
#define HANDLE_ERROR(ans) (handleError((ans), __FILE__, __LINE__))
inline void handleError(cudaError_t code, const char *file, int line)
{
   if (code != cudaSuccess) {
      fprintf(stderr,"CUDA assert: %s %s %d\n", cudaGetErrorString(code), file, line);
      exit(EXIT_FAILURE);
   }
}

#define HANDLE_KERNEL_ERROR(...) \
do {                                         \
    __VA_ARGS__;                             \
    HANDLE_ERROR( cudaPeekAtLastError() );   \
    HANDLE_ERROR( cudaDeviceSynchronize() ); \
} while(0)

#define RUN_KERNEL(kernel, th1, th2, ...) HANDLE_KERNEL_ERROR( kernel<<<dimBlock, dimGrid>>>(__VA_ARGS__) )
#endif

// Constants

#ifndef COMPUTE_ON_CPU
__constant__ 
#endif
lbm_consts d_consts;

const ssize_t V[][9] = {
    { 1, 1, 1, 0, 0, 0,-1,-1,-1 },
    { 1, 0,-1, 1, 0,-1, 1, 0,-1 }
};
const double T[] = { 1./36, 1./9, 1./36, 1./9, 4./9, 1./9, 1./36, 1./9, 1./36 };


/**
 * Setup: cylindrical obstacle and velocity inlet with perturbation
 * Creation of a mask with boolean values, defining the shape of the obstacle.
 */
static void initObstacles(lbm_vars* vars)
{
    for (int x = 0; x < NX; x++) {
        for (int y = 0; y < NY; y++) {
            vars->obstacles[x][y] = false;
        }
    }
}

/**
 * Initial velocity profile: almost zero, with a slight perturbation to trigger
 * the instability.
 */
static void initVelocity(lbm_vars* vars)
{
    for (int d = 0; d < 2; d++) {
        for (int x = 0; x < NX; x++) {
            for (int y = 0; y < NY; y++) {
                vars->vel[x][y][d] = 0; //(1-d) * ULB * (1 + 0.0001 * sin( y / (double)LY * 2 * M_PI) );
            }
        }
    }
}

static void initRho(lbm_vars* vars)
{
    for (int x = 0; x < NX; x++) {
        for (int y = 0; y < NY; y++) {
            vars->rho[x][y] = 1.0;
        }
    }
    vars->rho[NX/2][NY/2] = 2.0;
}

static void initCol(size_t* col, ssize_t v0)
{
    for (int f = 0, i = 0; f < 9 && i < 3; f++) {
        if (V[0][f] == v0) {
            col[i++] = f;
        }
    }
}

static void initOpp(size_t* opp)
{
    for (int f = 0; f < 9; f++) {
        for (int g = 0; g < 9; g++) {
            if (V[0][f] == -V[0][g] && V[1][f] == -V[1][g]) {
                opp[f] = g;
                break;
            }
        }
    }
}

#define EQUILIBRIUM_BODY(v, t) \
    do {                                                                 \
        double usqr = 3./2 * ( SQUARE(u[0]) + SQUARE(u[1]) );            \
                                                                         \
        for (int f = 0; f < 9; f++) {                                    \
            double cu = 3 * ( v[0][f] * u[0] + v[1][f] * u[1] );         \
            feq[f] = rho * t[f] * ( 1 + cu + 0.5 * SQUARE(cu) - usqr );  \
        }                                                                \
    } while(0)

#define GPU_EQUILIBRIUM_BODY(v, t) \
    do {                                                                                \
        double usqr = __dmul_rn(3./2, __dadd_rn( GPU_SQUARE(u[0]), GPU_SQUARE(u[1]) )); \
                                                                                        \
        for (int f = 0; f < 9; f++) {                                                   \
            double cu = 3 * ( v[0][f] * u[0] + v[1][f] * u[1] );                        \
            feq[f] = rho * t[f] * ( 1 + cu + 0.5 * SQUARE(cu) - usqr );                 \
        }                                                                               \
    } while(0)


#ifndef COMPUTE_ON_CPU
__host__ 
#endif
static void h_equilibrium(double* feq, double rho, double* u)
{
    EQUILIBRIUM_BODY(V, T);
}

#ifndef COMPUTE_ON_CPU
__device__ 
#endif
static void d_equilibrium(double* feq, double rho, double* u)
{
#ifdef COMPUTE_ON_CPU
    EQUILIBRIUM_BODY(d_consts.v, d_consts.t);
#else
    GPU_EQUILIBRIUM_BODY(d_consts.v, d_consts.t);
#endif
}

#ifndef COMPUTE_ON_CPU
__device__ 
#endif
static void macroscopic(double* fin, double* rho, double* u)
{
    
    *rho = u[0] = u[1] = 0;

    for (int f = 0; f < 9; f++) {
        *rho += fin[f];

        u[0] += d_consts.v[0][f] * fin[f];
        u[1] += d_consts.v[1][f] * fin[f];
    }
    
    u[0] /= *rho;
    u[1] /= *rho;
}

#ifdef COMPUTE_ON_CPU
void lbm_computation(lbm_vars *d_vars, int x, int y) 
{
#else
__global__ void lbm_computation(lbm_vars *d_vars)
{
    for (int idx = threadIdx.x + blockIdx.x * blockDim.x; idx < NX*NY; idx += blockDim.x * gridDim.x) {
        int x, y;
        INDEX_2D_FROM_1D(x, y, idx);
#endif

//    // Right wall: outflow condition.
//    if (x == NX-1) {
//        for (int i = 0; i < 3; i++) {
//            int f = d_consts.col[2][i];
//            d_vars->fin[NX-1][y][f] = d_vars->fin[NX-2][y][f];
//        }
//    }

    // Compute macroscopic variables, density and velocity
    macroscopic(d_vars->fin[x][y], &d_vars->rho[x][y], d_vars->u[x][y]);
    
//    if (x == 0) {
//        // Left wall: inflow condition
//        for (size_t d = 0; d < 2; d++) {
//            d_vars->u[0][y][d] = d_vars->vel[0][y][d];
//        }        

//        // Calculate the density
//        double s2 = 0, s3 = 0;
//        for (size_t i = 0; i < 3; i++) {
//            s2 += d_vars->fin[0][y][d_consts.col[1][i]];
//            s3 += d_vars->fin[0][y][d_consts.col[2][i]];
//        }
//        d_vars->rho[0][y] = 1./(1 - d_vars->u[0][y][0]) * (s2 + 2*s3);
//    }
   
    // Compute equilibrium
    d_equilibrium(d_vars->feq[x][y], d_vars->rho[x][y], d_vars->u[x][y]);

//    if (x == 0) {
//        for (size_t i = 0; i < 3; i++) {
//            size_t f = d_consts.col[0][i];
//            d_vars->fin[0][y][f] = d_vars->feq[0][y][f] + d_vars->fin[0][y][d_consts.opp[f]] - d_vars->feq[0][y][d_consts.opp[f]];
//        }
//    }
    
    for (size_t f = 0; f < 9; f++) {
        if (d_vars->obstacles[x][y]) {
            // Bounce-back condition for obstacle
            d_vars->fout[x][y][f] = d_vars->fin[x][y][d_consts.opp[f]];
        } else {
            // Collision step
#ifdef COMPUTE_ON_CPU
            d_vars->fout[x][y][f] = d_vars->fin[x][y][f] - OMEGA * (d_vars->fin[x][y][f] - d_vars->feq[x][y][f]);
#else
            d_vars->fout[x][y][f] = __dadd_rn(__dmul_rn(-OMEGA, __dadd_rn(d_vars->fin[x][y][f], - d_vars->feq[x][y][f])), d_vars->fin[x][y][f]);
#endif
        }
    }

#ifndef COMPUTE_ON_CPU
    }
#endif
}

#ifdef COMPUTE_ON_CPU
void lbm_streaming(lbm_vars *d_vars, int x, int y)
{
#else
__global__ void lbm_streaming(lbm_vars *d_vars)
{
    for (int idx = threadIdx.x + blockIdx.x * blockDim.x; idx < NX*NY; idx += blockDim.x * gridDim.x) {
        int x, y;
        INDEX_2D_FROM_1D(x, y, idx);
#endif
    // Streaming step
    for (size_t f = 0; f < 9; f++) {
        ssize_t x_dst = (x + d_consts.v[0][f] + NX) % NX ;
        ssize_t y_dst = (y + d_consts.v[1][f] + NY) % NY ;
        d_vars->fin[x_dst][y_dst][f] = d_vars->fout[x][y][f];
    }

#ifndef COMPUTE_ON_CPU
    }
#endif
}

void output_variables(char* filename, double var[NX][NY][9]) 
{
    FILE* file = fopen(filename, "w");
    for (size_t x = 0; x < NX; x++) {
        for (size_t y = 0; y < NY; y++) {
            for (size_t f = 0; f < 9; ++f) {
                fprintf(file, "%64.60f\n", var[x][y][f]);
            }
        }
    }
    fclose(file);
}

void output_image(char* filename, double u[NX][NY][2]) 
{
    pgm_image* pgm = pgm_create(NX, NY);

    for (size_t x = 0; x < NX; x++) {
        for (size_t y = 0; y < NY; y++) {
            double vel = sqrt( SQUARE(u[x][y][0]) + SQUARE(u[x][y][1]) );
            int color =  255 * min(vel * 10, 1.0);
{'fraction_non_alphanumeric': 0.17952342487883682, 'fraction_numerical': 0.0248889337641357, 'mean_word_length': 1.8708695652173912, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 97, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 46, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '9689', 'n_tokens_mistral': 8169, 'n_tokens_neox': 7789, 'n_words': 1836}