1 #ifndef LIBFLATARRAY_EXAMPLES_LBM_UPDATE_LBM_OBJECT_ORIENTED_H
2 #define LIBFLATARRAY_EXAMPLES_LBM_UPDATE_LBM_OBJECT_ORIENTED_H
3
4 /**
5 * Copyright 2013 Andreas Schäfer
6 *
7 * Distributed under the Boost Software License, Version 1.0. (See accompanying
8 * file LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
9 */
10
11 #include "util.h"
12
13 #define GET_COMP(X, Y, Z, DIR) \
14 gridOld[(Z) * dimX * dimY + (Y) * dimX + (X)].DIR
15
16 #define SET_COMP(DIR) \
17 gridNew[(z) * dimX * dimY + (y) * dimX + (x)].DIR
18
update_lbm_object_oriented(int dimX,int dimY,int dimZ,CellLBM * gridOld,CellLBM * gridNew)19 __global__ void update_lbm_object_oriented(int dimX, int dimY, int dimZ, CellLBM *gridOld, CellLBM *gridNew)
20 {
21 int x = blockIdx.x * blockDim.x + threadIdx.x + 2;
22 int y = blockIdx.y * blockDim.y + threadIdx.y + 2;
23 int z = 2;
24
25 #pragma unroll 10
26 for (; z < (dimZ - 2); z += 1) {
27
28 #define SQR(X) ((X)*(X))
29 const double omega = 1.0/1.7;
30 const double omega_trm = 1.0 - omega;
31 const double omega_w0 = 3.0 * 1.0 / 3.0 * omega;
32 const double omega_w1 = 3.0*1.0/18.0*omega;
33 const double omega_w2 = 3.0*1.0/36.0*omega;
34 const double one_third = 1.0 / 3.0;
35 double velX, velY, velZ;
36
37 velX =
38 GET_COMP(x-1,y,z,E) + GET_COMP(x-1,y-1,z,NE) +
39 GET_COMP(x-1,y+1,z,SE) + GET_COMP(x-1,y,z-1,TE) +
40 GET_COMP(x-1,y,z+1,BE);
41 velY = GET_COMP(x,y-1,z,N) + GET_COMP(x+1,y-1,z,NW) +
42 GET_COMP(x,y-1,z-1,TN) + GET_COMP(x,y-1,z+1,BN);
43 velZ = GET_COMP(x,y,z-1,T) + GET_COMP(x,y+1,z-1,TS) +
44 GET_COMP(x+1,y,z-1,TW);
45
46 const double rho =
47 GET_COMP(x,y,z,C) + GET_COMP(x,y+1,z,S) +
48 GET_COMP(x+1,y,z,W) + GET_COMP(x,y,z+1,B) +
49 GET_COMP(x+1,y+1,z,SW) + GET_COMP(x,y+1,z+1,BS) +
50 GET_COMP(x+1,y,z+1,BW) + velX + velY + velZ;
51 velX = velX
52 - GET_COMP(x+1,y,z,W) - GET_COMP(x+1,y-1,z,NW)
53 - GET_COMP(x+1,y+1,z,SW) - GET_COMP(x+1,y,z-1,TW)
54 - GET_COMP(x+1,y,z+1,BW);
55 velY = velY
56 + GET_COMP(x-1,y-1,z,NE) - GET_COMP(x,y+1,z,S)
57 - GET_COMP(x+1,y+1,z,SW) - GET_COMP(x-1,y+1,z,SE)
58 - GET_COMP(x,y+1,z-1,TS) - GET_COMP(x,y+1,z+1,BS);
59 velZ = velZ+GET_COMP(x,y-1,z-1,TN) + GET_COMP(x-1,y,z-1,TE) - GET_COMP(x,y,z+1,B) - GET_COMP(x,y-1,z+1,BN) - GET_COMP(x,y+1,z+1,BS) - GET_COMP(x+1,y,z+1,BW) - GET_COMP(x-1,y,z+1,BE);
60
61 // density = rho;
62 // velocityX = velX;
63 // velocityY = velY;
64 // velocityZ = velZ;
65
66 const double dir_indep_trm = one_third*rho - 0.5*( velX*velX + velY*velY + velZ*velZ );
67
68 SET_COMP(C)=omega_trm * GET_COMP(x,y,z,C) + omega_w0*( dir_indep_trm );
69
70 SET_COMP(NW)=omega_trm * GET_COMP(x+1,y-1,z,NW) +
71 omega_w2*( dir_indep_trm - ( velX-velY ) + 1.5*SQR( velX-velY ) );
72 SET_COMP(SE)=omega_trm * GET_COMP(x-1,y+1,z,SE) +
73 omega_w2*( dir_indep_trm + ( velX-velY ) + 1.5*SQR( velX-velY ) );
74 SET_COMP(NE)=omega_trm * GET_COMP(x-1,y-1,z,NE) +
75 omega_w2*( dir_indep_trm + ( velX+velY ) + 1.5*SQR( velX+velY ) );
76 SET_COMP(SW)=omega_trm * GET_COMP(x+1,y+1,z,SW) +
77 omega_w2*( dir_indep_trm - ( velX+velY ) + 1.5*SQR( velX+velY ) );
78
79 SET_COMP(TW)=omega_trm * GET_COMP(x+1,y,z-1,TW) + omega_w2*( dir_indep_trm - ( velX-velZ ) + 1.5*SQR( velX-velZ ) );
80 SET_COMP(BE)=omega_trm * GET_COMP(x-1,y,z+1,BE) + omega_w2*( dir_indep_trm + ( velX-velZ ) + 1.5*SQR( velX-velZ ) );
81 SET_COMP(TE)=omega_trm * GET_COMP(x-1,y,z-1,TE) + omega_w2*( dir_indep_trm + ( velX+velZ ) + 1.5*SQR( velX+velZ ) );
82 SET_COMP(BW)=omega_trm * GET_COMP(x+1,y,z+1,BW) + omega_w2*( dir_indep_trm - ( velX+velZ ) + 1.5*SQR( velX+velZ ) );
83
84 SET_COMP(TS)=omega_trm * GET_COMP(x,y+1,z-1,TS) + omega_w2*( dir_indep_trm - ( velY-velZ ) + 1.5*SQR( velY-velZ ) );
85 SET_COMP(BN)=omega_trm * GET_COMP(x,y-1,z+1,BN) + omega_w2*( dir_indep_trm + ( velY-velZ ) + 1.5*SQR( velY-velZ ) );
86 SET_COMP(TN)=omega_trm * GET_COMP(x,y-1,z-1,TN) + omega_w2*( dir_indep_trm + ( velY+velZ ) + 1.5*SQR( velY+velZ ) );
87 SET_COMP(BS)=omega_trm * GET_COMP(x,y+1,z+1,BS) + omega_w2*( dir_indep_trm - ( velY+velZ ) + 1.5*SQR( velY+velZ ) );
88
89 SET_COMP(N)=omega_trm * GET_COMP(x,y-1,z,N) + omega_w1*( dir_indep_trm + velY + 1.5*SQR(velY));
90 SET_COMP(S)=omega_trm * GET_COMP(x,y+1,z,S) + omega_w1*( dir_indep_trm - velY + 1.5*SQR(velY));
91 SET_COMP(E)=omega_trm * GET_COMP(x-1,y,z,E) + omega_w1*( dir_indep_trm + velX + 1.5*SQR(velX));
92 SET_COMP(W)=omega_trm * GET_COMP(x+1,y,z,W) + omega_w1*( dir_indep_trm - velX + 1.5*SQR(velX));
93 SET_COMP(T)=omega_trm * GET_COMP(x,y,z-1,T) + omega_w1*( dir_indep_trm + velZ + 1.5*SQR(velZ));
94 SET_COMP(B)=omega_trm * GET_COMP(x,y,z+1,B) + omega_w1*( dir_indep_trm - velZ + 1.5*SQR(velZ));
95 }
96 }
97
98 #undef GET_COMP
99 #undef SET_COMP
100
101 class benchmark_lbm_cuda_object_oriented : public benchmark_lbm_cuda_basic
102 {
103 public:
name()104 virtual std::string name()
105 {
106 return "lbm_cuda_object_oriented";
107 }
108
109 protected:
update(dim3 dimGrid,dim3 dimBlock,int dimX,int dimY,int dimZ,double * devGridOld,double * devGridNew)110 void update(dim3 dimGrid, dim3 dimBlock, int dimX, int dimY, int dimZ, double *devGridOld, double *devGridNew)
111 {
112 update_lbm_object_oriented<<<dimGrid, dimBlock>>>(
113 dimX, dimY, dimZ,
114 reinterpret_cast<CellLBM*>(devGridOld),
115 reinterpret_cast<CellLBM*>(devGridNew));
116 }
117 };
118
119 #endif
120