mirror of
https://github.com/smyalygames/FiniteVolumeGPU.git
synced 2025-05-18 06:24:13 +02:00
Refactoring
This commit is contained in:
André R. Brodtkorb
parent
d9eb72d78c
commit
0671bd747a
@ -28,27 +28,18 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
__device__
|
||||
void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float F[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
const float g_, const float dx_, const float dt_) {
|
||||
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
//Compute fluxes along the x axis
|
||||
{
|
||||
int j=ty;
|
||||
const int l = j + 1; //Skip ghost cells
|
||||
for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
const int k = i;
|
||||
|
||||
for (int j=threadIdx.y; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
|
||||
for (int i=threadIdx.x; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
// Q at interface from the right and left
|
||||
const float3 Qp = make_float3(Q[0][l][k+1],
|
||||
Q[1][l][k+1],
|
||||
Q[2][l][k+1]);
|
||||
const float3 Qm = make_float3(Q[0][l][k],
|
||||
Q[1][l][k],
|
||||
Q[2][l][k]);
|
||||
const float3 Qp = make_float3(Q[0][j][i+1],
|
||||
Q[1][j][i+1],
|
||||
Q[2][j][i+1]);
|
||||
const float3 Qm = make_float3(Q[0][j][i],
|
||||
Q[1][j][i],
|
||||
Q[2][j][i]);
|
||||
|
||||
// Computed flux
|
||||
const float3 flux = FORCE_1D_flux(Qm, Qp, g_, dx_, dt_);
|
||||
@ -65,27 +56,19 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
*/
|
||||
__device__
|
||||
void computeFluxG(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float G[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
const float g_, const float dy_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
//Compute fluxes along the y axis
|
||||
for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
const int l = j;
|
||||
{
|
||||
int i=tx;
|
||||
const int k = i + 1; //Skip ghost cells
|
||||
|
||||
for (int j=threadIdx.y; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
for (int i=threadIdx.x; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
|
||||
// Q at interface from the right and left
|
||||
// Note that we swap hu and hv
|
||||
const float3 Qp = make_float3(Q[0][l+1][k],
|
||||
Q[2][l+1][k],
|
||||
Q[1][l+1][k]);
|
||||
const float3 Qm = make_float3(Q[0][l][k],
|
||||
Q[2][l][k],
|
||||
Q[1][l][k]);
|
||||
const float3 Qp = make_float3(Q[0][j+1][i],
|
||||
Q[2][j+1][i],
|
||||
Q[1][j+1][i]);
|
||||
const float3 Qm = make_float3(Q[0][j][i],
|
||||
Q[2][j][i],
|
||||
Q[1][j][i]);
|
||||
|
||||
// Computed flux
|
||||
// Note that we swap back
|
||||
@ -120,7 +103,7 @@ __global__ void FORCEKernel(
|
||||
const unsigned int vars = 3;
|
||||
|
||||
__shared__ float Q[3][h+2][w+2];
|
||||
__shared__ float F[3][h+1][w+1];
|
||||
__shared__ float F[3][h+2][w+2];
|
||||
|
||||
//Read into shared memory
|
||||
readBlock<w, h, gc>( h0_ptr_, h0_pitch_, Q[0], nx_+2, ny_+2);
|
||||
|
||||
@ -31,24 +31,20 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
__device__
|
||||
void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float F[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
const float g_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
{
|
||||
const int j=ty;
|
||||
const int l = j + 1; //Skip ghost cells
|
||||
for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
const int k = i;
|
||||
|
||||
const float3 Q_l = make_float3(Q[0][l][k ], Q[1][l][k ], Q[2][l][k ]);
|
||||
const float3 Q_r = make_float3(Q[0][l][k+1], Q[1][l][k+1], Q[2][l][k+1]);
|
||||
|
||||
for (int j=threadIdx.y; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
|
||||
for (int i=threadIdx.x; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
// Q at interface from the right and left
|
||||
const float3 Q_r = make_float3(Q[0][j][i+1],
|
||||
Q[1][j][i+1],
|
||||
Q[2][j][i+1]);
|
||||
const float3 Q_l = make_float3(Q[0][j][i],
|
||||
Q[1][j][i],
|
||||
Q[2][j][i]);
|
||||
|
||||
// Computed flux
|
||||
const float3 flux = HLL_flux(Q_l, Q_r, g_);
|
||||
|
||||
//Write to shared memory
|
||||
F[0][j][i] = flux.x;
|
||||
F[1][j][i] = flux.y;
|
||||
F[2][j][i] = flux.z;
|
||||
@ -65,27 +61,23 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
*/
|
||||
__device__
|
||||
void computeFluxG(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float G[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
const float g_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
const int l = j;
|
||||
{
|
||||
const int i=tx;
|
||||
const int k = i + 1; //Skip ghost cells
|
||||
|
||||
//NOte that hu and hv are swapped ("transposing" the domain)!
|
||||
const float3 Q_l = make_float3(Q[0][l ][k], Q[2][l ][k], Q[1][l ][k]);
|
||||
const float3 Q_r = make_float3(Q[0][l+1][k], Q[2][l+1][k], Q[1][l+1][k]);
|
||||
//Compute fluxes along the y axis
|
||||
for (int j=threadIdx.y; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
for (int i=threadIdx.x; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
|
||||
// Q at interface from the right and left
|
||||
// Note that we swap hu and hv
|
||||
const float3 Q_r = make_float3(Q[0][j+1][i],
|
||||
Q[2][j+1][i],
|
||||
Q[1][j+1][i]);
|
||||
const float3 Q_l = make_float3(Q[0][j][i],
|
||||
Q[2][j][i],
|
||||
Q[1][j][i]);
|
||||
|
||||
// Computed flux
|
||||
const float3 flux = HLL_flux(Q_l, Q_r, g_);
|
||||
|
||||
//Write to shared memory
|
||||
//Note that we here swap hu and hv back to the original
|
||||
const float3 flux = HLL_flux(Q_l, Q_r, g_);
|
||||
G[0][j][i] = flux.x;
|
||||
G[1][j][i] = flux.z;
|
||||
G[2][j][i] = flux.y;
|
||||
@ -128,7 +120,7 @@ __global__ void HLLKernel(
|
||||
|
||||
//Shared memory variables
|
||||
__shared__ float Q[3][h+2][w+2];
|
||||
__shared__ float F[3][h+1][w+1];
|
||||
__shared__ float F[3][h+2][w+2];
|
||||
|
||||
//Read into shared memory
|
||||
readBlock<w, h, gc>( h0_ptr_, h0_pitch_, Q[0], nx_+2, ny_+2);
|
||||
|
||||
@ -33,33 +33,26 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
__device__
|
||||
void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float Qx[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float g_, const float dx_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
{
|
||||
const int j=ty;
|
||||
const int l = j + 2; //Skip ghost cells
|
||||
for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
const int k = i + 1;
|
||||
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 float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
const float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float3 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
|
||||
@ -85,34 +78,27 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
*/
|
||||
__device__
|
||||
void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float Qy[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float g_, const float dy_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
const int l = j + 1;
|
||||
{
|
||||
int i=tx;
|
||||
const int k = i + 2; //Skip ghost cells
|
||||
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 float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
const float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float3 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
|
||||
@ -163,8 +149,8 @@ __global__ void HLL2Kernel(
|
||||
|
||||
//Shared memory variables
|
||||
__shared__ float Q[3][h+4][w+4];
|
||||
__shared__ float Qx[3][h+2][w+2];
|
||||
__shared__ float F[3][h+1][w+1];
|
||||
__shared__ float Qx[3][h+4][w+4];
|
||||
__shared__ float F[3][h+4][w+4];
|
||||
|
||||
//Read into shared memory
|
||||
readBlock<w, h, gc>( h0_ptr_, h0_pitch_, Q[0], nx_+2, ny_+2);
|
||||
|
||||
@ -96,6 +96,39 @@ void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
|
||||
|
||||
|
||||
__device__ void minmodSlopeX(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
const float theta_) {
|
||||
//Reconstruct slopes along x axis
|
||||
for (int p=0; p<3; ++p) {
|
||||
{
|
||||
const int j = threadIdx.y+2;
|
||||
for (int i=threadIdx.x+1; i<BLOCK_WIDTH+3; i+=BLOCK_WIDTH) {
|
||||
Qx[p][j-2][i-1] = minmodSlope(Q[p][j][i-1], Q[p][j][i], Q[p][j][i+1], theta_);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Reconstructs a minmod slope for a whole block along the ordinate
|
||||
*/
|
||||
__device__ void minmodSlopeY(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
const float theta_) {
|
||||
//Reconstruct slopes along y axis
|
||||
for (int p=0; p<3; ++p) {
|
||||
const int i = threadIdx.x + 2;
|
||||
for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+3; j+=BLOCK_HEIGHT) {
|
||||
{
|
||||
Qy[p][j-1][i-2] = minmodSlope(Q[p][j-1][i], Q[p][j][i], Q[p][j+1][i], theta_);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
|
||||
*/
|
||||
@ -159,8 +192,8 @@ __global__ void KP07Kernel(
|
||||
|
||||
|
||||
//Reconstruct slopes along x and axis
|
||||
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_);
|
||||
minmodSlopeY<w, h, gc, vars>(Q, Qy, theta_);
|
||||
minmodSlopeX(Q, Qx, theta_);
|
||||
minmodSlopeY(Q, Qy, theta_);
|
||||
__syncthreads();
|
||||
|
||||
|
||||
|
||||
@ -31,77 +31,65 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
__device__
|
||||
void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float Qx[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float g_, const float dx_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
int j=ty;
|
||||
const int l = j + 2; //Skip ghost cells
|
||||
for (int i=tx; 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 float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float3 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
|
||||
const float3 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, g_) - F_func(Q_lr, g_));
|
||||
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 float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
|
||||
// Compute flux based on prediction
|
||||
const float3 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, g_);
|
||||
|
||||
//Write to shared memory
|
||||
F[0][j][i] = flux.x;
|
||||
F[1][j][i] = flux.y;
|
||||
F[2][j][i] = flux.z;
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float3 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
|
||||
const float3 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, g_) - F_func(Q_lr, g_));
|
||||
|
||||
// Compute flux based on prediction
|
||||
const float3 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, g_);
|
||||
|
||||
//Write to shared memory
|
||||
F[0][j][i] = flux.x;
|
||||
F[1][j][i] = flux.y;
|
||||
F[2][j][i] = flux.z;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
__device__
|
||||
void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
|
||||
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float Qy[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float g_, const float dy_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
const int l = j + 1;
|
||||
{
|
||||
int i=tx;
|
||||
const int k = i + 2; //Skip ghost cells
|
||||
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 float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
const float3 Q_rl = make_float3(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]);
|
||||
const float3 Q_rr = make_float3(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]);
|
||||
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
const float3 Q_ll = make_float3(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]);
|
||||
const float3 Q_lr = make_float3(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]);
|
||||
|
||||
//Evolve half a timestep (predictor step)
|
||||
const float3 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
|
||||
@ -153,8 +141,8 @@ __global__ void KP07DimsplitKernel(
|
||||
|
||||
//Shared memory variables
|
||||
__shared__ float Q[3][h+4][w+4];
|
||||
__shared__ float Qx[3][h+2][w+2];
|
||||
__shared__ float F[3][h+1][w+1];
|
||||
__shared__ float Qx[3][h+4][w+4];
|
||||
__shared__ float F[3][h+4][w+4];
|
||||
|
||||
|
||||
|
||||
|
||||
@ -34,23 +34,15 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
__device__
|
||||
void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float g_, const float dx_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
{
|
||||
int j=ty;
|
||||
const int l = j + 2; //Skip ghost cells
|
||||
for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
|
||||
const int k = i + 1;
|
||||
|
||||
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) {
|
||||
// Q at interface from the right and left
|
||||
const float3 Ql2 = make_float3(Q[0][l][k-1], Q[1][l][k-1], Q[2][l][k-1]);
|
||||
const float3 Ql1 = make_float3(Q[0][l][k ], Q[1][l][k ], Q[2][l][k ]);
|
||||
const float3 Qr1 = make_float3(Q[0][l][k+1], Q[1][l][k+1], Q[2][l][k+1]);
|
||||
const float3 Qr2 = make_float3(Q[0][l][k+2], Q[1][l][k+2], Q[2][l][k+2]);
|
||||
const float3 Ql2 = make_float3(Q[0][j][i-1], Q[1][j][i-1], Q[2][j][i-1]);
|
||||
const float3 Ql1 = make_float3(Q[0][j][i ], Q[1][j][i ], Q[2][j][i ]);
|
||||
const float3 Qr1 = make_float3(Q[0][j][i+1], Q[1][j][i+1], Q[2][j][i+1]);
|
||||
const float3 Qr2 = make_float3(Q[0][j][i+2], Q[1][j][i+2], Q[2][j][i+2]);
|
||||
|
||||
// Computed flux
|
||||
const float3 flux = WAF_1D_flux(Ql2, Ql1, Qr1, Qr2, g_, dx_, dt_);
|
||||
@ -73,24 +65,16 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
*/
|
||||
__device__
|
||||
void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
|
||||
float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
|
||||
const float g_, const float dy_, const float dt_) {
|
||||
//Index of thread within block
|
||||
const int tx = threadIdx.x;
|
||||
const int ty = threadIdx.y;
|
||||
|
||||
//Compute fluxes along the y axis
|
||||
for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
|
||||
const int l = j + 1;
|
||||
{
|
||||
int i=tx;
|
||||
const int k = i + 2; //Skip ghost cells
|
||||
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) {
|
||||
// Q at interface from the right and left
|
||||
// Note that we swap hu and hv
|
||||
const float3 Ql2 = make_float3(Q[0][l-1][k], Q[2][l-1][k], Q[1][l-1][k]);
|
||||
const float3 Ql1 = make_float3(Q[0][l ][k], Q[2][l ][k], Q[1][l ][k]);
|
||||
const float3 Qr1 = make_float3(Q[0][l+1][k], Q[2][l+1][k], Q[1][l+1][k]);
|
||||
const float3 Qr2 = make_float3(Q[0][l+2][k], Q[2][l+2][k], Q[1][l+2][k]);
|
||||
const float3 Ql2 = make_float3(Q[0][j-1][i], Q[2][j-1][i], Q[1][j-1][i]);
|
||||
const float3 Ql1 = make_float3(Q[0][j ][i], Q[2][j ][i], Q[1][j ][i]);
|
||||
const float3 Qr1 = make_float3(Q[0][j+1][i], Q[2][j+1][i], Q[1][j+1][i]);
|
||||
const float3 Qr2 = make_float3(Q[0][j+2][i], Q[2][j+2][i], Q[1][j+2][i]);
|
||||
|
||||
// Computed flux
|
||||
// Note that we swap back
|
||||
@ -138,7 +122,7 @@ __global__ void WAFKernel(
|
||||
|
||||
//Shared memory variables
|
||||
__shared__ float Q[3][h+4][w+4];
|
||||
__shared__ float F[3][h+1][w+1];
|
||||
__shared__ float F[3][h+4][w+4];
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user