mirror of
https://github.com/smyalygames/FiniteVolumeGPU.git
synced 2025-11-29 17:28:03 +01:00
Refactoring
This commit is contained in:
André R. Brodtkorb
@@ -27,97 +27,95 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
__device__
|
||||
void computeFluxF(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qx[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float F[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
const float gamma_, const float dx_, const float dt_) {
|
||||
int j=threadIdx.y;
|
||||
const int l = j + 2; //Skip ghost cells
|
||||
for (int i=threadIdx.x; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
const int k = i + 1;
|
||||
// Reconstruct point values of Q at the left and right hand side
|
||||
// of the cell for both the left (i) and right (i+1) cell
|
||||
const float4 Q_rl = make_float4(Q[0][l][k+1] - 0.5f*Qx[0][j][i+1],
|
||||
Q[1][l][k+1] - 0.5f*Qx[1][j][i+1],
|
||||
Q[2][l][k+1] - 0.5f*Qx[2][j][i+1],
|
||||
Q[3][l][k+1] - 0.5f*Qx[3][j][i+1]);
|
||||
const float4 Q_rr = make_float4(Q[0][l][k+1] + 0.5f*Qx[0][j][i+1],
|
||||
Q[1][l][k+1] + 0.5f*Qx[1][j][i+1],
|
||||
Q[2][l][k+1] + 0.5f*Qx[2][j][i+1],
|
||||
Q[3][l][k+1] + 0.5f*Qx[3][j][i+1]);
|
||||
|
||||
const float4 Q_ll = make_float4(Q[0][l][k] - 0.5f*Qx[0][j][i],
|
||||
Q[1][l][k] - 0.5f*Qx[1][j][i],
|
||||
Q[2][l][k] - 0.5f*Qx[2][j][i],
|
||||
Q[3][l][k] - 0.5f*Qx[3][j][i]);
|
||||
const float4 Q_lr = make_float4(Q[0][l][k] + 0.5f*Qx[0][j][i],
|
||||
Q[1][l][k] + 0.5f*Qx[1][j][i],
|
||||
Q[2][l][k] + 0.5f*Qx[2][j][i],
|
||||
Q[3][l][k] + 0.5f*Qx[3][j][i]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float4 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
|
||||
const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
|
||||
float Qx[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float F[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float gamma_, const float dx_, const float dt_) {
|
||||
for (int j=threadIdx.y; j<BLOCK_HEIGHT+4; j+=BLOCK_HEIGHT) {
|
||||
for (int i=threadIdx.x+1; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
|
||||
// Reconstruct point values of Q at the left and right hand side
|
||||
// of the cell for both the left (i) and right (i+1) cell
|
||||
const float4 Q_rl = make_float4(Q[0][j][i+1] - 0.5f*Qx[0][j][i+1],
|
||||
Q[1][j][i+1] - 0.5f*Qx[1][j][i+1],
|
||||
Q[2][j][i+1] - 0.5f*Qx[2][j][i+1],
|
||||
Q[3][j][i+1] - 0.5f*Qx[3][j][i+1]);
|
||||
const float4 Q_rr = make_float4(Q[0][j][i+1] + 0.5f*Qx[0][j][i+1],
|
||||
Q[1][j][i+1] + 0.5f*Qx[1][j][i+1],
|
||||
Q[2][j][i+1] + 0.5f*Qx[2][j][i+1],
|
||||
Q[3][j][i+1] + 0.5f*Qx[3][j][i+1]);
|
||||
|
||||
// Compute flux based on prediction
|
||||
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
|
||||
|
||||
//Write to shared memory
|
||||
F[0][j][i] = flux.x;
|
||||
F[1][j][i] = flux.y;
|
||||
F[2][j][i] = flux.z;
|
||||
F[3][j][i] = flux.w;
|
||||
const float4 Q_ll = make_float4(Q[0][j][i] - 0.5f*Qx[0][j][i],
|
||||
Q[1][j][i] - 0.5f*Qx[1][j][i],
|
||||
Q[2][j][i] - 0.5f*Qx[2][j][i],
|
||||
Q[3][j][i] - 0.5f*Qx[3][j][i]);
|
||||
const float4 Q_lr = make_float4(Q[0][j][i] + 0.5f*Qx[0][j][i],
|
||||
Q[1][j][i] + 0.5f*Qx[1][j][i],
|
||||
Q[2][j][i] + 0.5f*Qx[2][j][i],
|
||||
Q[3][j][i] + 0.5f*Qx[3][j][i]);
|
||||
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float4 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
|
||||
const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
|
||||
|
||||
// Compute flux based on prediction
|
||||
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
|
||||
|
||||
//Write to shared memory
|
||||
F[0][j][i] = flux.x;
|
||||
F[1][j][i] = flux.y;
|
||||
F[2][j][i] = flux.z;
|
||||
F[3][j][i] = flux.w;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
__device__
|
||||
void computeFluxG(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qy[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float G[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
const float gamma_, const float dy_, const float dt_) {
|
||||
int i=threadIdx.x;
|
||||
const int k = i + 2; //Skip ghost cells
|
||||
for (int j=threadIdx.y; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
const int l = j + 1;
|
||||
// Reconstruct point values of Q at the left and right hand side
|
||||
// of the cell for both the left (i) and right (i+1) cell
|
||||
//NOte that hu and hv are swapped ("transposing" the domain)!
|
||||
const float4 Q_rl = make_float4(Q[0][l+1][k] - 0.5f*Qy[0][j+1][i],
|
||||
Q[2][l+1][k] - 0.5f*Qy[2][j+1][i],
|
||||
Q[1][l+1][k] - 0.5f*Qy[1][j+1][i],
|
||||
Q[3][l+1][k] - 0.5f*Qy[3][j+1][i]);
|
||||
const float4 Q_rr = make_float4(Q[0][l+1][k] + 0.5f*Qy[0][j+1][i],
|
||||
Q[2][l+1][k] + 0.5f*Qy[2][j+1][i],
|
||||
Q[1][l+1][k] + 0.5f*Qy[1][j+1][i],
|
||||
Q[3][l+1][k] + 0.5f*Qy[3][j+1][i]);
|
||||
float Qy[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float G[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float gamma_, const float dy_, const float dt_) {
|
||||
for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
|
||||
for (int i=threadIdx.x; i<BLOCK_WIDTH+4; i+=BLOCK_WIDTH) {
|
||||
// Reconstruct point values of Q at the left and right hand side
|
||||
// of the cell for both the left (i) and right (i+1) cell
|
||||
//NOte that hu and hv are swapped ("transposing" the domain)!
|
||||
const float4 Q_rl = make_float4(Q[0][j+1][i] - 0.5f*Qy[0][j+1][i],
|
||||
Q[2][j+1][i] - 0.5f*Qy[2][j+1][i],
|
||||
Q[1][j+1][i] - 0.5f*Qy[1][j+1][i],
|
||||
Q[3][j+1][i] - 0.5f*Qy[3][j+1][i]);
|
||||
const float4 Q_rr = make_float4(Q[0][j+1][i] + 0.5f*Qy[0][j+1][i],
|
||||
Q[2][j+1][i] + 0.5f*Qy[2][j+1][i],
|
||||
Q[1][j+1][i] + 0.5f*Qy[1][j+1][i],
|
||||
Q[3][j+1][i] + 0.5f*Qy[3][j+1][i]);
|
||||
|
||||
const float4 Q_ll = make_float4(Q[0][l][k] - 0.5f*Qy[0][j][i],
|
||||
Q[2][l][k] - 0.5f*Qy[2][j][i],
|
||||
Q[1][l][k] - 0.5f*Qy[1][j][i],
|
||||
Q[3][l][k] - 0.5f*Qy[3][j][i]);
|
||||
const float4 Q_lr = make_float4(Q[0][l][k] + 0.5f*Qy[0][j][i],
|
||||
Q[2][l][k] + 0.5f*Qy[2][j][i],
|
||||
Q[1][l][k] + 0.5f*Qy[1][j][i],
|
||||
Q[3][l][k] + 0.5f*Qy[3][j][i]);
|
||||
const float4 Q_ll = make_float4(Q[0][j][i] - 0.5f*Qy[0][j][i],
|
||||
Q[2][j][i] - 0.5f*Qy[2][j][i],
|
||||
Q[1][j][i] - 0.5f*Qy[1][j][i],
|
||||
Q[3][j][i] - 0.5f*Qy[3][j][i]);
|
||||
const float4 Q_lr = make_float4(Q[0][j][i] + 0.5f*Qy[0][j][i],
|
||||
Q[2][j][i] + 0.5f*Qy[2][j][i],
|
||||
Q[1][j][i] + 0.5f*Qy[1][j][i],
|
||||
Q[3][j][i] + 0.5f*Qy[3][j][i]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float4 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
|
||||
const float4 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
|
||||
|
||||
// Compute flux based on prediction
|
||||
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
|
||||
|
||||
//Write to shared memory
|
||||
//Note that we here swap hu and hv back to the original
|
||||
G[0][j][i] = flux.x;
|
||||
G[1][j][i] = flux.z;
|
||||
G[2][j][i] = flux.y;
|
||||
G[3][j][i] = flux.w;
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float4 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
|
||||
const float4 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
|
||||
|
||||
// Compute flux based on prediction
|
||||
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
|
||||
|
||||
//Write to shared memory
|
||||
//Note that we here swap hu and hv back to the original
|
||||
G[0][j][i] = flux.x;
|
||||
G[1][j][i] = flux.z;
|
||||
G[2][j][i] = flux.y;
|
||||
G[3][j][i] = flux.w;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
|
||||
*/
|
||||
@@ -150,8 +148,8 @@ __global__ void KP07DimsplitKernel(
|
||||
|
||||
//Shared memory variables
|
||||
__shared__ float Q[4][h+4][w+4];
|
||||
__shared__ float Qx[4][h+2][w+2];
|
||||
__shared__ float F[4][h+1][w+1];
|
||||
__shared__ float Qx[4][h+4][w+4];
|
||||
__shared__ float F[4][h+4][w+4];
|
||||
|
||||
|
||||
|
||||
@@ -162,14 +160,12 @@ __global__ void KP07DimsplitKernel(
|
||||
readBlock<w, h, gc>( E0_ptr_, E0_pitch_, Q[3], nx_+4, ny_+4);
|
||||
__syncthreads();
|
||||
|
||||
|
||||
//Fix boundary conditions
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[0], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, -1, 1>(Q[1], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[1], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, -1>(Q[2], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[3], nx_, ny_);
|
||||
__syncthreads();
|
||||
|
||||
|
||||
|
||||
//Step 0 => evolve x first, then y
|
||||
@@ -186,12 +182,11 @@ __global__ void KP07DimsplitKernel(
|
||||
|
||||
//Set boundary conditions
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[0], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, -1, 1>(Q[1], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[1], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, -1>(Q[2], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[3], nx_, ny_);
|
||||
__syncthreads();
|
||||
|
||||
|
||||
//Compute fluxes along the y axis and evolve
|
||||
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_);
|
||||
__syncthreads();
|
||||
@@ -200,7 +195,7 @@ __global__ void KP07DimsplitKernel(
|
||||
__syncthreads();
|
||||
|
||||
evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
|
||||
__syncthreads();
|
||||
__syncthreads();
|
||||
|
||||
}
|
||||
//Step 1 => evolve y first, then x
|
||||
@@ -208,14 +203,16 @@ __global__ void KP07DimsplitKernel(
|
||||
//Compute fluxes along the y axis and evolve
|
||||
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_);
|
||||
__syncthreads();
|
||||
|
||||
computeFluxG(Q, Qx, F, gamma_, dy_, dt_);
|
||||
__syncthreads();
|
||||
|
||||
evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
|
||||
__syncthreads();
|
||||
|
||||
|
||||
//Set boundary conditions
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[0], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, -1, 1>(Q[1], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[1], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, -1>(Q[2], nx_, ny_);
|
||||
noFlowBoundary<w, h, gc, 1, 1>(Q[3], nx_, ny_);
|
||||
__syncthreads();
|
||||
@@ -223,8 +220,10 @@ __global__ void KP07DimsplitKernel(
|
||||
//Compute fluxes along the x axis and evolve
|
||||
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_);
|
||||
__syncthreads();
|
||||
|
||||
computeFluxF(Q, Qx, F, gamma_, dx_, dt_);
|
||||
__syncthreads();
|
||||
|
||||
evolveF<w, h, gc, vars>(Q, F, dx_, dt_);
|
||||
__syncthreads();
|
||||
}
|
||||
|
||||
@@ -151,86 +151,91 @@ inline __device__ void writeBlock(float* ptr_, int pitch_,
|
||||
|
||||
template<int block_width, int block_height, int ghost_cells, int scale_east_west=1, int scale_north_south=1>
|
||||
__device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) {
|
||||
const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells;
|
||||
const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells;
|
||||
|
||||
const int i = threadIdx.x + ghost_cells;
|
||||
const int j = threadIdx.y + ghost_cells;
|
||||
|
||||
|
||||
// West boundary
|
||||
if (ti == ghost_cells) {
|
||||
Q[j][i-1] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && ti == ghost_cells + 1) {
|
||||
Q[j][i-3] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && ti == ghost_cells + 2) {
|
||||
Q[j][i-5] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 4 && ti == ghost_cells + 3) {
|
||||
Q[j][i-7] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 5 && ti == ghost_cells + 4) {
|
||||
Q[j][i-9] = scale_east_west*Q[j][i];
|
||||
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+= block_height) {
|
||||
const int i = threadIdx.x + ghost_cells;
|
||||
const int ti = blockDim.x*blockIdx.x + i;
|
||||
const int tj = blockDim.y*blockIdx.y + j;
|
||||
|
||||
// West boundary
|
||||
if (ti == ghost_cells) {
|
||||
Q[j][i-1] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && ti == ghost_cells + 1) {
|
||||
Q[j][i-3] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && ti == ghost_cells + 2) {
|
||||
Q[j][i-5] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 4 && ti == ghost_cells + 3) {
|
||||
Q[j][i-7] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 5 && ti == ghost_cells + 4) {
|
||||
Q[j][i-9] = scale_east_west*Q[j][i];
|
||||
}
|
||||
|
||||
|
||||
|
||||
// East boundary
|
||||
if (ti == nx_ + ghost_cells - 1) {
|
||||
Q[j][i+1] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && ti == nx_ + ghost_cells - 2) {
|
||||
Q[j][i+3] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 3) {
|
||||
Q[j][i+5] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 4 && ti == nx_ + ghost_cells - 4) {
|
||||
Q[j][i+7] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 5 && ti == nx_ + ghost_cells - 5) {
|
||||
Q[j][i+9] = scale_east_west*Q[j][i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// East boundary
|
||||
if (ti == nx_ + ghost_cells - 1) {
|
||||
Q[j][i+1] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && ti == nx_ + ghost_cells - 2) {
|
||||
Q[j][i+3] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 3) {
|
||||
Q[j][i+5] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 4) {
|
||||
Q[j][i+7] = scale_east_west*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 5) {
|
||||
Q[j][i+9] = scale_east_west*Q[j][i];
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
// South boundary
|
||||
if (tj == ghost_cells) {
|
||||
Q[j-1][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && tj == ghost_cells + 1) {
|
||||
Q[j-3][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && tj == ghost_cells + 2) {
|
||||
Q[j-5][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 4 && tj == ghost_cells + 3) {
|
||||
Q[j-7][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 5 && tj == ghost_cells + 4) {
|
||||
Q[j-9][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
|
||||
|
||||
|
||||
// North boundary
|
||||
if (tj == ny_ + ghost_cells - 1) {
|
||||
Q[j+1][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && tj == ny_ + ghost_cells - 2) {
|
||||
Q[j+3][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 3) {
|
||||
Q[j+5][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 4) {
|
||||
Q[j+7][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 5) {
|
||||
Q[j+9][i] = scale_north_south*Q[j][i];
|
||||
|
||||
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+= block_width) {
|
||||
const int j = threadIdx.y + ghost_cells;
|
||||
const int ti = blockDim.x*blockIdx.x + i;
|
||||
const int tj = blockDim.y*blockIdx.y + j;
|
||||
|
||||
// South boundary
|
||||
if (tj == ghost_cells) {
|
||||
Q[j-1][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && tj == ghost_cells + 1) {
|
||||
Q[j-3][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && tj == ghost_cells + 2) {
|
||||
Q[j-5][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 4 && tj == ghost_cells + 3) {
|
||||
Q[j-7][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 5 && tj == ghost_cells + 4) {
|
||||
Q[j-9][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
|
||||
|
||||
|
||||
// North boundary
|
||||
if (tj == ny_ + ghost_cells - 1) {
|
||||
Q[j+1][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 2 && tj == ny_ + ghost_cells - 2) {
|
||||
Q[j+3][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 3) {
|
||||
Q[j+5][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 4 && tj == ny_ + ghost_cells - 4) {
|
||||
Q[j+7][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
if (ghost_cells >= 5 && tj == ny_ + ghost_cells - 5) {
|
||||
Q[j+9][i] = scale_north_south*Q[j][i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -247,21 +252,14 @@ __device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2
|
||||
|
||||
template<int block_width, int block_height, int ghost_cells, int vars>
|
||||
__device__ void evolveF(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
float F[vars][block_height+1][block_width+1],
|
||||
float F[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
const float dx_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
const int i = tx + ghost_cells; //Skip local ghost cells
|
||||
const int j = ty + ghost_cells;
|
||||
|
||||
//Index of cell within domain
|
||||
//const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells; //Skip global ghost cells, i.e., +1
|
||||
//const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells;
|
||||
//if (ti > ghost_cells-1 && ti < nx_+ghost_cells && tj > ghost_cells-1 && tj < ny_+ghost_cells) {
|
||||
for (int var=0; var < vars; ++var) {
|
||||
Q[var][j][i] = Q[var][j][i] + (F[var][ty][tx] - F[var][ty][tx+1]) * dt_ / dx_;
|
||||
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) {
|
||||
for (int i=threadIdx.x+1; i<block_width+2*ghost_cells; i+=block_width) {
|
||||
Q[var][j][i] = Q[var][j][i] + (F[var][j][i-1] - F[var][j][i]) * dt_ / dx_;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -275,17 +273,14 @@ __device__ void evolveF(float Q[vars][block_height+2*ghost_cells][block_width+2*
|
||||
*/
|
||||
template<int block_width, int block_height, int ghost_cells, int vars>
|
||||
__device__ void evolveG(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
float G[vars][block_height+1][block_width+1],
|
||||
float G[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
const float dy_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
const int i = tx + ghost_cells; //Skip local ghost cells, i.e., +1
|
||||
const int j = ty + ghost_cells;
|
||||
|
||||
for (int var=0; var < vars; ++var) {
|
||||
Q[var][j][i] = Q[var][j][i] + (G[var][ty][tx] - G[var][ty+1][tx]) * dt_ / dy_;
|
||||
for (int j=threadIdx.y+1; j<block_height+2*ghost_cells; j+=block_height) {
|
||||
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) {
|
||||
Q[var][j][i] = Q[var][j][i] + (G[var][j-1][i] - G[var][j][i]) * dt_ / dy_;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -328,5 +323,6 @@ __device__ void memset(float Q[vars][shmem_height][shmem_width], float value) {
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -47,21 +47,15 @@ __device__ __inline__ float minmodSlope(float left, float center, float right, f
|
||||
* Reconstructs a minmod slope for a whole block along the abscissa
|
||||
*/
|
||||
template<int block_width, int block_height, int ghost_cells, int vars>
|
||||
__device__ void minmodSlopeX(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
float Qx[vars][block_height+2*(ghost_cells-1)][block_width+2*(ghost_cells-1)],
|
||||
__device__ void minmodSlopeX(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
float Qx[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
const float theta_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
const int j = ty;
|
||||
const int l = j + ghost_cells; //Skip ghost cells
|
||||
|
||||
//Reconstruct slopes along x axis
|
||||
for (int i=tx; i<block_width+2*(ghost_cells-1); i+=block_width) {
|
||||
const int k = i + 1;
|
||||
for (int p=0; p<vars; ++p) {
|
||||
Qx[p][j][i] = minmodSlope(Q[p][l][k-1], Q[p][l][k], Q[p][l][k+1], theta_);
|
||||
for (int p=0; p<vars; ++p) {
|
||||
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) {
|
||||
for (int i=threadIdx.x+1; i<block_width+3; i+=block_width) {
|
||||
Qx[p][j][i] = minmodSlope(Q[p][j][i-1], Q[p][j][i], Q[p][j][i+1], theta_);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -71,21 +65,15 @@ __device__ void minmodSlopeX(float Q[vars][block_height+2*ghost_cells][block_wi
|
||||
* Reconstructs a minmod slope for a whole block along the ordinate
|
||||
*/
|
||||
template<int block_width, int block_height, int ghost_cells, int vars>
|
||||
__device__ void minmodSlopeY(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
float Qy[vars][block_height+2*(ghost_cells-1)][block_width+2*(ghost_cells-1)],
|
||||
__device__ void minmodSlopeY(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
float Qy[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
|
||||
const float theta_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
const int i = tx;
|
||||
const int k = i + ghost_cells; //Skip ghost cells
|
||||
|
||||
//Reconstruct slopes along y axis
|
||||
for (int j=ty; j<block_height+2*(ghost_cells-1); j+=block_height) {
|
||||
const int l = j + 1;
|
||||
for (int p=0; p<vars; ++p) {
|
||||
Qy[p][j][i] = minmodSlope(Q[p][l-1][k], Q[p][l][k], Q[p][l+1][k], theta_);
|
||||
for (int p=0; p<vars; ++p) {
|
||||
for (int j=threadIdx.y+1; j<block_height+3; j+=block_height) {
|
||||
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) {
|
||||
Qy[p][j][i] = minmodSlope(Q[p][j-1][i], Q[p][j][i], Q[p][j+1][i], theta_);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user