smyalygames/FiniteVolumeGPU
1
0
Fork 0
mirror of https://github.com/smyalygames/FiniteVolumeGPU.git synced 2025-11-29 17:28:03 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity

Refactoring

Browse Source
This commit is contained in:
André R. Brodtkorb
2018-09-04 16:33:50 +02:00
parent 58f281d724
commit 10d8e26108
30 changed files with 774 additions and 840 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
GPUSimulators/cuda/EulerCommon.h Normal file
View File

@@ -0,0 +1,59 @@
/*
This OpenCL kernel implements the Kurganov-Petrova numerical scheme
for the shallow water equations, described in
A. Kurganov & Guergana Petrova
A Second-Order Well-Balanced Positivity Preserving Central-Upwind
Scheme for the Saint-Venant System Communications in Mathematical
Sciences, 5 (2007), 133-160.
Copyright (C) 2016 SINTEF ICT
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/>.
*/
#pragma once
__device__ float pressure(float4 Q, float gamma) {
const float rho = Q.x;
const float rho_u = Q.y;
const float rho_v = Q.z;
const float E = Q.w;
return (gamma-1.0f)*(E-0.5f*(rho_u*rho_u + rho_v*rho_v)/rho);
}
__device__ float4 F_func(const float4 Q, float gamma) {
const float rho = Q.x;
const float rho_u = Q.y;
const float rho_v = Q.z;
const float E = Q.w;
const float u = rho_u/rho;
const float P = pressure(Q, gamma);
float4 F;
F.x = rho_u;
F.y = rho_u*u + P;
F.z = rho_v*u;
F.w = u*(E+P);
return F;
}

474
GPUSimulators/cuda/SWECommon.h Normal file
View File

@@ -0,0 +1,474 @@
/*
This OpenCL kernel implements the Kurganov-Petrova numerical scheme
for the shallow water equations, described in
A. Kurganov & Guergana Petrova
A Second-Order Well-Balanced Positivity Preserving Central-Upwind
Scheme for the Saint-Venant System Communications in Mathematical
Sciences, 5 (2007), 133-160.
Copyright (C) 2016 SINTEF ICT
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/>.
*/
#pragma once
#include "limiters.h"
__device__ float3 F_func(const float3 Q, const float g) {
float3 F;
F.x = Q.y; //hu
F.y = Q.y*Q.y / Q.x + 0.5f*g*Q.x*Q.x; //hu*hu/h + 0.5f*g*h*h;
F.z = Q.y*Q.z / Q.x; //hu*hv/h;
return F;
}
/**
* Superbee flux limiter for WAF.
* Related to superbee limiter so that WAF_superbee(r, c) = 1 - (1-|c|)*superbee(r)
* @param r_ the ratio of upwind change (see Toro 2001, p. 203/204)
* @param c_ the courant number for wave k, dt*S_k/dx
*/
__device__ float WAF_superbee(float r_, float c_) {
// r <= 0.0
if (r_ <= 0.0f) {
return 1.0f;
}
// 0.0 <= r <= 1/2
else if (r_ <= 0.5f) {
return 1.0f - 2.0f*(1.0f - fabsf(c_))*r_;
}
// 1/2 <= r <= 1
else if (r_ <= 1.0f) {
return fabs(c_);
}
// 1 <= r <= 2
else if (r_ <= 2.0f) {
return 1.0f - (1.0f - fabsf(c_))*r_;
}
// r >= 2
else {
return 2.0f*fabsf(c_) - 1.0f;
}
}
__device__ float WAF_albada(float r_, float c_) {
if (r_ <= 0.0f) {
return 1.0f;
}
else {
return 1.0f - (1.0f - fabsf(c_)) * r_ * (1.0f + r_) / (1.0f + r_*r_);
}
}
__device__ float WAF_minbee(float r_, float c_) {
r_ = fmaxf(-1.0f, fminf(2.0f, r_));
if (r_ <= 0.0f) {
return 1.0f;
}
if (r_ >= 0.0f && r_ <= 1.0f) {
return 1.0f - (1.0f - fabsf(c_)) * r_;
}
else {
return fabsf(c_);
}
}
__device__ float WAF_minmod(float r_, float c_) {
return 1.0f - (1.0f - fabsf(c_)) * fmaxf(0.0f, fminf(1.0f, r_));
}
__device__ float limiterToWAFLimiter(float r_, float c_) {
return 1.0f - (1.0f - fabsf(c_))*r_;
}
// Compute h in the "star region", h^dagger
__device__ __inline__ float computeHStar(float h_l, float h_r, float u_l, float u_r, float c_l, float c_r, float g_) {
//This estimate for the h* gives rise to spurious oscillations.
//return 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
const float h_tmp = 0.5f * (c_l + c_r) + 0.25f * (u_l - u_r);
return h_tmp*h_tmp / g_;
}
/**
* Weighted average flux (Toro 2001, p 200) for interface {i+1/2}
* @param r_ The flux limiter parameter (see Toro 2001, p. 203)
* @param Q_l2 Q_{i-1}
* @param Q_l1 Q_{i}
* @param Q_r1 Q_{i+1}
* @param Q_r2 Q_{i+2}
*/
__device__ float3 WAF_1D_flux(const float3 Q_l2, const float3 Q_l1, const float3 Q_r1, const float3 Q_r2, const float g_, const float dx_, const float dt_) {
const float h_l = Q_l1.x;
const float h_r = Q_r1.x;
const float h_l2 = Q_l2.x;
const float h_r2 = Q_r2.x;
// Calculate velocities
const float u_l = Q_l1.y / h_l;
const float u_r = Q_r1.y / h_r;
const float u_l2 = Q_l2.y / h_l2;
const float u_r2 = Q_r2.y / h_r2;
const float v_l = Q_l1.z / h_l;
const float v_r = Q_r1.z / h_r;
const float v_l2 = Q_l2.z / h_l2;
const float v_r2 = Q_r2.z / h_r2;
// Estimate the potential wave speeds
const float c_l = sqrt(g_*h_l);
const float c_r = sqrt(g_*h_r);
const float c_l2 = sqrt(g_*h_l2);
const float c_r2 = sqrt(g_*h_r2);
// Compute h in the "star region", h^dagger
const float h_dag_l = computeHStar(h_l2, h_l, u_l2, u_l, c_l2, c_l, g_);
const float h_dag = computeHStar( h_l, h_r, u_l, u_r, c_l, c_r, g_);
const float h_dag_r = computeHStar( h_r, h_r2, u_r, u_r2, c_r, c_r2, g_);
const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag ) ) / h_l;
const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag ) ) / h_r;
const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
// Compute wave speed estimates
const float S_l = u_l - c_l*q_l;
const float S_r = u_r + c_r*q_r;
const float S_star = ( S_l*h_r*(u_r - S_r) - S_r*h_l*(u_l - S_l) ) / ( h_r*(u_r - S_r) - h_l*(u_l - S_l) );
const float3 Q_star_l = h_l * (S_l - u_l) / (S_l - S_star) * make_float3(1.0, S_star, v_l);
const float3 Q_star_r = h_r * (S_r - u_r) / (S_r - S_star) * make_float3(1.0, S_star, v_r);
// Estimate the fluxes in the four regions
const float3 F_1 = F_func(Q_l1, g_);
const float3 F_4 = F_func(Q_r1, g_);
const float3 F_2 = F_1 + S_l*(Q_star_l - Q_l1);
const float3 F_3 = F_4 + S_r*(Q_star_r - Q_r1);
//const float3 F_2 = F_func(Q_star_l, g_);
//const float3 F_3 = F_func(Q_star_r, g_);
// Compute the courant numbers for the waves
const float c_1 = S_l * dt_ / dx_;
const float c_2 = S_star * dt_ / dx_;
const float c_3 = S_r * dt_ / dx_;
// Compute the "upwind change" vectors for the i-3/2 and i+3/2 interfaces
const float eps = 1.0e-6f;
const float r_1 = desingularize( (c_1 > 0.0f) ? (h_dag_l - h_l2) : (h_dag_r - h_r), eps) / desingularize((h_dag - h_l), eps);
const float r_2 = desingularize( (c_2 > 0.0f) ? (v_l - v_l2) : (v_r2 - v_r), eps ) / desingularize((v_r - v_l), eps);
const float r_3 = desingularize( (c_3 > 0.0f) ? (h_l - h_dag_l) : (h_r2 - h_dag_r), eps ) / desingularize((h_r - h_dag), eps);
// Compute the limiter
// We use h for the nonlinear waves, and v for the middle shear wave
const float A_1 = copysign(1.0f, c_1) * limiterToWAFLimiter(generalized_minmod(r_1, 1.9f), c_1);
const float A_2 = copysign(1.0f, c_2) * limiterToWAFLimiter(generalized_minmod(r_2, 1.9f), c_2);
const float A_3 = copysign(1.0f, c_3) * limiterToWAFLimiter(generalized_minmod(r_3, 1.9f), c_3);
//Average the fluxes
const float3 flux = 0.5f*( F_1 + F_4 )
- 0.5f*( A_1 * (F_2 - F_1)
+ A_2 * (F_3 - F_2)
+ A_3 * (F_4 - F_3) );
return flux;
}
/**
* Central upwind flux function
*/
__device__ float3 CentralUpwindFlux(const float3 Qm, float3 Qp, const float g) {
const float3 Fp = F_func(Qp, g);
const float up = Qp.y / Qp.x; // hu / h
const float cp = sqrt(g*Qp.x); // sqrt(g*h)
const float3 Fm = F_func(Qm, g);
const float um = Qm.y / Qm.x; // hu / h
const float cm = sqrt(g*Qm.x); // sqrt(g*h)
const float am = min(min(um-cm, up-cp), 0.0f); // largest negative wave speed
const float ap = max(max(um+cm, up+cp), 0.0f); // largest positive wave speed
return ((ap*Fm - am*Fp) + ap*am*(Qp-Qm))/(ap-am);
}
/**
* Godunovs centered scheme (Toro 2001, p 165)
*/
__device__ float3 GodC_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
const float3 F_l = F_func(Q_l, g_);
const float3 F_r = F_func(Q_r, g_);
const float3 Q_godc = 0.5f*(Q_l + Q_r) + (dt_/dx_)*(F_l - F_r);
return F_func(Q_godc, g_);
}
/**
* Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 180)
*/
__device__ float3 HLL_flux(const float3 Q_l, const float3 Q_r, const float g_) {
const float h_l = Q_l.x;
const float h_r = Q_r.x;
// Calculate velocities
const float u_l = Q_l.y / h_l;
const float u_r = Q_r.y / h_r;
// Estimate the potential wave speeds
const float c_l = sqrt(g_*h_l);
const float c_r = sqrt(g_*h_r);
// Compute h in the "star region", h^dagger
const float h_dag = 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag / (h_l*h_l) ) );
const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag / (h_r*h_r) ) );
const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
// Compute wave speed estimates
const float S_l = u_l - c_l*q_l;
const float S_r = u_r + c_r*q_r;
//Upwind selection
if (S_l >= 0.0f) {
return F_func(Q_l, g_);
}
else if (S_r <= 0.0f) {
return F_func(Q_r, g_);
}
//Or estimate flux in the star region
else {
const float3 F_l = F_func(Q_l, g_);
const float3 F_r = F_func(Q_r, g_);
const float3 flux = (S_r*F_l - S_l*F_r + S_r*S_l*(Q_r - Q_l)) / (S_r-S_l);
return flux;
}
}
/**
* Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 181)
*/
__device__ float3 HLLC_flux(const float3 Q_l, const float3 Q_r, const float g_) {
const float h_l = Q_l.x;
const float h_r = Q_r.x;
// Calculate velocities
const float u_l = Q_l.y / h_l;
const float u_r = Q_r.y / h_r;
// Estimate the potential wave speeds
const float c_l = sqrt(g_*h_l);
const float c_r = sqrt(g_*h_r);
// Compute h in the "star region", h^dagger
const float h_dag = 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag / (h_l*h_l) ) );
const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag / (h_r*h_r) ) );
const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
// Compute wave speed estimates
const float S_l = u_l - c_l*q_l;
const float S_r = u_r + c_r*q_r;
const float S_star = ( S_l*h_r*(u_r - S_r) - S_r*h_l*(u_l - S_l) ) / ( h_r*(u_r - S_r) - h_l*(u_l - S_l) );
const float3 F_l = F_func(Q_l, g_);
const float3 F_r = F_func(Q_r, g_);
//Upwind selection
if (S_l >= 0.0f) {
return F_l;
}
else if (S_r <= 0.0f) {
return F_r;
}
//Or estimate flux in the "left star" region
else if (S_l <= 0.0f && 0.0f <=S_star) {
const float v_l = Q_l.z / h_l;
const float3 Q_star_l = h_l * (S_l - u_l) / (S_l - S_star) * make_float3(1, S_star, v_l);
const float3 flux = F_l + S_l*(Q_star_l - Q_l);
return flux;
}
//Or estimate flux in the "righ star" region
else if (S_star <= 0.0f && 0.0f <=S_r) {
const float v_r = Q_r.z / h_r;
const float3 Q_star_r = h_r * (S_r - u_r) / (S_r - S_star) * make_float3(1, S_star, v_r);
const float3 flux = F_r + S_r*(Q_star_r - Q_r);
return flux;
}
else {
return make_float3(-99999.9f, -99999.9f, -99999.9f); //Something wrong here
}
}
/**
* Lax-Friedrichs flux (Toro 2001, p 163)
*/
__device__ float3 LxF_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
const float3 F_l = F_func(Q_l, g_);
const float3 F_r = F_func(Q_r, g_);
return 0.5f*(F_l + F_r) + (dx_/(2.0f*dt_))*(Q_l - Q_r);
}
/**
* Lax-Friedrichs extended to 2D
*/
__device__ float3 LxF_2D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
const float3 F_l = F_func(Q_l, g_);
const float3 F_r = F_func(Q_r, g_);
//Note numerical diffusion for 2D here (0.25)
return 0.5f*(F_l + F_r) + (dx_/(4.0f*dt_))*(Q_l - Q_r);
}
/**
* Richtmeyer / Two-step Lax-Wendroff flux (Toro 2001, p 164)
*/
__device__ float3 LxW2_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
const float3 F_l = F_func(Q_l, g_);
const float3 F_r = F_func(Q_r, g_);
const float3 Q_lw2 = 0.5f*(Q_l + Q_r) + (dt_/(2.0f*dx_))*(F_l - F_r);
return F_func(Q_lw2, g_);
}
/**
* First Ordered Centered (Toro 2001, p.163)
*/
__device__ float3 FORCE_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
const float3 F_lf = LxF_1D_flux(Q_l, Q_r, g_, dx_, dt_);
const float3 F_lw2 = LxW2_1D_flux(Q_l, Q_r, g_, dx_, dt_);
return 0.5f*(F_lf + F_lw2);
}

153
GPUSimulators/cuda/SWE_FORCE.cu Normal file
View File

@@ -0,0 +1,153 @@
/*
This OpenCL kernel implements the classical Lax-Friedrichs scheme
for the shallow water equations, with edge fluxes.
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
/**
* Computes the flux along the x axis for all faces
*/
__device__
void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
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;
// 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]);
// Computed flux
const float3 flux = FORCE_1D_flux(Qm, Qp, g_, dx_, dt_);
F[0][j][i] = flux.x;
F[1][j][i] = flux.y;
F[2][j][i] = flux.z;
}
}
}
/**
* Computes the flux along the y axis for all faces
*/
__device__
void computeFluxG(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
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
// 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]);
// Computed flux
// Note that we swap back
const float3 flux = FORCE_1D_flux(Qm, Qp, g_, dy_, dt_);
G[0][j][i] = flux.x;
G[1][j][i] = flux.z;
G[2][j][i] = flux.y;
}
}
}
extern "C" {
__global__ void FORCEKernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
__shared__ float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
//Read into shared memory
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+2, BLOCK_HEIGHT+2, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+2, ny_+2);
__syncthreads();
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
__syncthreads();
//Compute flux along x, and evolve
computeFluxF(Q, F, g_, dx_, dt_);
__syncthreads();
evolveF1(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
__syncthreads();
//Compute flux along y, and evolve
computeFluxG(Q, F, g_, dy_, dt_);
__syncthreads();
evolveG1(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
//Write to main memory
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_);
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_);
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_);
}
} // extern "C"

161
GPUSimulators/cuda/SWE_HLL.cu Normal file
View File

@@ -0,0 +1,161 @@
/*
This GPU kernel implements the HLL flux
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
/**
* Computes the flux along the x axis for all faces
*/
__device__
void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
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]);
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;
}
}
}
/**
* Computes the flux along the y axis for all faces
*/
__device__
void computeFluxG(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
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]);
// 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
G[0][j][i] = flux.x;
G[1][j][i] = flux.z;
G[2][j][i] = flux.y;
}
}
}
extern "C" {
__global__ void HLLKernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
//Shared memory variables
__shared__ float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
//Read into shared memory
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+2, BLOCK_HEIGHT+2, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+2, ny_+2);
__syncthreads();
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
__syncthreads();
//Compute F flux
computeFluxF(Q, F, g_);
__syncthreads();
evolveF1(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
__syncthreads();
//Compute G flux
computeFluxG(Q, F, g_);
__syncthreads();
evolveG1(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
// Write to main memory for all internal cells
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_);
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_);
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_);
}
} // extern "C"

232
GPUSimulators/cuda/SWE_HLL2.cu Normal file
View File

@@ -0,0 +1,232 @@
/*
This OpenCL kernel implements the second order HLL flux
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
#include "limiters.h"
/**
* Computes the flux along the x axis for all faces
*/
__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],
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;
// 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_));
// Compute flux based on prediction
const float3 flux = HLL_flux(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;
}
}
}
/**
* Computes the flux along the x axis for all faces
*/
__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],
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
// 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_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]);
//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_));
const float3 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (F_func(Q_ll, g_) - F_func(Q_lr, g_));
// Compute flux based on prediction
const float3 flux = HLL_flux(Q_l_bar, Q_r_bar, g_);
//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;
}
}
}
extern "C" {
__global__ void HLL2Kernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
float theta_,
int step_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
//Shared memory variables
__shared__ float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4];
__shared__ float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
//Read into shared memory
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+4, BLOCK_HEIGHT+4, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+4, ny_+4);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Step 0 => evolve x first, then y
if (step_ == 0) {
//Compute fluxes along the x axis and evolve
minmodSlopeX(Q, Qx, theta_);
__syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the y axis and evolve
minmodSlopeY(Q, Qx, theta_);
__syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
}
//Step 1 => evolve y first, then x
else {
//Compute fluxes along the y axis and evolve
minmodSlopeY(Q, Qx, theta_);
__syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the x axis and evolve
minmodSlopeX(Q, Qx, theta_);
__syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
}
// Write to main memory for all internal cells
writeBlock2(h1_ptr_, h1_pitch_,
hu1_ptr_, hu1_pitch_,
hv1_ptr_, hv1_pitch_,
Q, nx_, ny_);
}
} // extern "C"

209
GPUSimulators/cuda/SWE_KP07.cu Normal file
View File

@@ -0,0 +1,209 @@
/*
This OpenCL kernel implements the Kurganov-Petrova numerical scheme
for the shallow water equations, described in
A. Kurganov & Guergana Petrova
A Second-Order Well-Balanced Positivity Preserving Central-Upwind
Scheme for the Saint-Venant System Communications in Mathematical
Sciences, 5 (2007), 133-160.
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
#include "limiters.h"
__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],
const float g_) {
//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;
// Q at interface from the right and left
const float3 Qp = 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 Qm = 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 ]);
// Computed flux
const float3 flux = CentralUpwindFlux(Qm, Qp, g_);
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],
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 + 1;
{
int i=tx;
const int k = i + 2; //Skip ghost cells
// 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] - 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 Qm = 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]);
// Computed flux
// Note that we swap back
const float3 flux = CentralUpwindFlux(Qm, Qp, g_);
G[0][j][i] = flux.x;
G[1][j][i] = flux.z;
G[2][j][i] = flux.y;
}
}
}
/**
* This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
*/
extern "C" {
__global__ void KP07Kernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
float theta_,
int step_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2;
//Shared memory variables
__shared__ float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4];
//The following slightly wastes memory, but enables us to reuse the
//funcitons in common.opencl
__shared__ float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
__shared__ float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
//Read into shared memory
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+4, BLOCK_HEIGHT+4, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+4, ny_+4);
__syncthreads();
//Fix boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Reconstruct slopes along x and axis
minmodSlopeX(Q, Qx, theta_);
minmodSlopeY(Q, Qy, theta_);
__syncthreads();
//Compute fluxes along the x and y axis
computeFluxF(Q, Qx, F, g_);
computeFluxG(Q, Qy, G, g_);
__syncthreads();
//Sum fluxes and advance in time for all internal cells
if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) {
const int i = tx + 2; //Skip local ghost cells, i.e., +2
const int j = ty + 2;
const float h1 = Q[0][j][i] + (F[0][ty][tx] - F[0][ty ][tx+1]) * dt_ / dx_
+ (G[0][ty][tx] - G[0][ty+1][tx ]) * dt_ / dy_;
const float hu1 = Q[1][j][i] + (F[1][ty][tx] - F[1][ty ][tx+1]) * dt_ / dx_
+ (G[1][ty][tx] - G[1][ty+1][tx ]) * dt_ / dy_;
const float hv1 = Q[2][j][i] + (F[2][ty][tx] - F[2][ty ][tx+1]) * dt_ / dx_
+ (G[2][ty][tx] - G[2][ty+1][tx ]) * dt_ / dy_;
float* const h_row = (float*) ((char*) h1_ptr_ + h1_pitch_*tj);
float* const hu_row = (float*) ((char*) hu1_ptr_ + hu1_pitch_*tj);
float* const hv_row = (float*) ((char*) hv1_ptr_ + hv1_pitch_*tj);
if (step_ == 0) {
//First step of RK2 ODE integrator
h_row[ti] = h1;
hu_row[ti] = hu1;
hv_row[ti] = hv1;
}
else if (step_ == 1) {
//Second step of RK2 ODE integrator
//First read Q^n
const float h_a = h_row[ti];
const float hu_a = hu_row[ti];
const float hv_a = hv_row[ti];
//Compute Q^n+1
const float h_b = 0.5f*(h_a + h1);
const float hu_b = 0.5f*(hu_a + hu1);
const float hv_b = 0.5f*(hv_a + hv1);
//Write to main memory
h_row[ti] = h_b;
hu_row[ti] = hu_b;
hv_row[ti] = hv_b;
}
}
}
} //extern "C"

224
GPUSimulators/cuda/SWE_KP07_dimsplit.cu Normal file
View File

@@ -0,0 +1,224 @@
/*
This OpenCL kernel implements the Kurganov-Petrova numerical scheme
for the shallow water equations, described in
A. Kurganov & Guergana Petrova
A Second-Order Well-Balanced Positivity Preserving Central-Upwind
Scheme for the Saint-Venant System Communications in Mathematical
Sciences, 5 (2007), 133-160.
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
#include "limiters.h"
__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],
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_));
// 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],
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
// 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_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]);
//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_));
const float3 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (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
//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;
}
}
}
/**
* This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
*/
extern "C" {
__global__ void KP07DimsplitKernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
float theta_,
int step_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
//Shared memory variables
__shared__ float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4];
__shared__ float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
//Read into shared memory
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+4, BLOCK_HEIGHT+4, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+4, ny_+4);
__syncthreads();
//Fix boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Step 0 => evolve x first, then y
if (step_ == 0) {
//Compute fluxes along the x axis and evolve
minmodSlopeX(Q, Qx, theta_);
__syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the y axis and evolve
minmodSlopeY(Q, Qx, theta_);
__syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
}
//Step 1 => evolve y first, then x
else {
//Compute fluxes along the y axis and evolve
minmodSlopeY(Q, Qx, theta_);
__syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the x axis and evolve
minmodSlopeX(Q, Qx, theta_);
__syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
}
// Write to main memory for all internal cells
writeBlock2(h1_ptr_, h1_pitch_,
hu1_ptr_, hu1_pitch_,
hv1_ptr_, hv1_pitch_,
Q, nx_, ny_);
}
} // extern "C"

158
GPUSimulators/cuda/SWE_LxF.cu Normal file
View File

@@ -0,0 +1,158 @@
/*
This OpenCL kernel implements the classical Lax-Friedrichs scheme
for the shallow water equations, with edge fluxes.
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
/**
* Computes the flux along the x axis for all faces
*/
template <int block_width, int block_height>
__device__
void computeFluxF(float Q[3][block_height+2][block_width+2],
float F[3][block_height][block_width+1],
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 + 1; //Skip ghost cells
for (int i=tx; i<block_width+1; i+=block_width) {
const int k = i;
// 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]);
// Computed flux
const float3 flux = LxF_2D_flux(Qm, Qp, g_, dx_, dt_);
F[0][j][i] = flux.x;
F[1][j][i] = flux.y;
F[2][j][i] = flux.z;
}
}
}
/**
* Computes the flux along the y axis for all faces
*/
template <int block_width, int block_height>
__device__
void computeFluxG(float Q[3][block_height+2][block_width+2],
float G[3][block_height+1][block_width],
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;
{
const int i=tx;
const int k = i + 1; //Skip ghost cells
// 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]);
// Computed flux
// Note that we swap back
const float3 flux = LxF_2D_flux(Qm, Qp, g_, dy_, dt_);
G[0][j][i] = flux.x;
G[1][j][i] = flux.z;
G[2][j][i] = flux.y;
}
}
}
extern "C" {
__global__
void LxFKernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
const int tx = threadIdx.x;
const int ty = threadIdx.y;
__shared__ float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2];
__shared__ float F[3][BLOCK_HEIGHT][BLOCK_WIDTH+1];
__shared__ float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH];
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+2, BLOCK_HEIGHT+2, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+2, ny_+2);
__syncthreads();
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the x and y axis
computeFluxF<BLOCK_WIDTH, BLOCK_HEIGHT>(Q, F, g_, dx_, dt_);
computeFluxG<BLOCK_WIDTH, BLOCK_HEIGHT>(Q, G, g_, dy_, dt_);
__syncthreads();
//Evolve for all cells
const int i = tx + 1; //Skip local ghost cells, i.e., +1
const int j = ty + 1;
Q[0][j][i] += (F[0][ty][tx] - F[0][ty ][tx+1]) * dt_ / dx_
+ (G[0][ty][tx] - G[0][ty+1][tx ]) * dt_ / dy_;
Q[1][j][i] += (F[1][ty][tx] - F[1][ty ][tx+1]) * dt_ / dx_
+ (G[1][ty][tx] - G[1][ty+1][tx ]) * dt_ / dy_;
Q[2][j][i] += (F[2][ty][tx] - F[2][ty ][tx+1]) * dt_ / dx_
+ (G[2][ty][tx] - G[2][ty+1][tx ]) * dt_ / dy_;
//Write to main memory
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_);
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_);
writeBlock<BLOCK_WIDTH+2, BLOCK_HEIGHT+2, 1, 1>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_);
}
} // extern "C"

198
GPUSimulators/cuda/SWE_WAF.cu Normal file
View File

@@ -0,0 +1,198 @@
/*
This OpenCL kernel implements the Kurganov-Petrova numerical scheme
for the shallow water equations, described in
A. Kurganov & Guergana Petrova
A Second-Order Well-Balanced Positivity Preserving Central-Upwind
Scheme for the Saint-Venant System Communications in Mathematical
Sciences, 5 (2007), 133-160.
Copyright (C) 2016 SINTEF ICT
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 "common.h"
#include "SWECommon.h"
/**
* Computes the flux along the x axis for all faces
*/
__device__
void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
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;
// 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]);
// Computed flux
const float3 flux = WAF_1D_flux(Ql2, Ql1, Qr1, Qr2, g_, dx_, dt_);
F[0][j][i] = flux.x;
F[1][j][i] = flux.y;
F[2][j][i] = flux.z;
}
}
}
/**
* Computes the flux along the y axis for all faces
*/
__device__
void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
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
// 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]);
// Computed flux
// Note that we swap back
const float3 flux = WAF_1D_flux(Ql2, Ql1, Qr1, Qr2, g_, dy_, dt_);
G[0][j][i] = flux.x;
G[1][j][i] = flux.z;
G[2][j][i] = flux.y;
}
}
}
extern "C" {
__global__ void WAFKernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_, int step_,
//Input h^n
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
//Shared memory variables
__shared__ float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4];
__shared__ float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1];
//Read into shared memory Q from global memory
float* Q_ptr[3] = {h0_ptr_, hu0_ptr_, hv0_ptr_};
int Q_pitch[3] = {h0_pitch_, hu0_pitch_, hv0_pitch_};
readBlock<3, BLOCK_WIDTH+4, BLOCK_HEIGHT+4, BLOCK_WIDTH, BLOCK_HEIGHT>(Q_ptr, Q_pitch, Q, nx_+4, ny_+4);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Step 0 => evolve x first, then y
if (step_ == 0) {
//Compute fluxes along the x axis and evolve
computeFluxF(Q, F, g_, dx_, dt_);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
//Fix boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the y axis and evolve
computeFluxG(Q, F, g_, dy_, dt_);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
}
//Step 1 => evolve y first, then x
else {
//Compute fluxes along the y axis and evolve
computeFluxG(Q, F, g_, dy_, dt_);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
__syncthreads();
//Fix boundary conditions
noFlowBoundary2(Q, nx_, ny_);
__syncthreads();
//Compute fluxes along the x axis and evolve
computeFluxF(Q, F, g_, dx_, dt_);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
__syncthreads();
}
// Write to main memory for all internal cells
writeBlock2(h1_ptr_, h1_pitch_,
hu1_ptr_, hu1_pitch_,
hv1_ptr_, hv1_pitch_,
Q, nx_, ny_);
}
} // extern "C"

394
GPUSimulators/cuda/common.h Normal file
View File

@@ -0,0 +1,394 @@
/*
This OpenCL kernel implements the Kurganov-Petrova numerical scheme
for the shallow water equations, described in
A. Kurganov & Guergana Petrova
A Second-Order Well-Balanced Positivity Preserving Central-Upwind
Scheme for the Saint-Venant System Communications in Mathematical
Sciences, 5 (2007), 133-160.
Copyright (C) 2016 SINTEF ICT
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/>.
*/
#pragma once
/**
* Float3 operators
*/
inline __device__ float3 operator*(const float a, const float3 b) {
return make_float3(a*b.x, a*b.y, a*b.z);
}
inline __device__ float3 operator/(const float3 a, const float b) {
return make_float3(a.x/b, a.y/b, a.z/b);
}
inline __device__ float3 operator-(const float3 a, const float3 b) {
return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
}
inline __device__ float3 operator+(const float3 a, const float3 b) {
return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
}
inline __device__ __host__ float clamp(const float f, const float a, const float b) {
return fmaxf(a, fminf(f, b));
}
__device__ float desingularize(float x_, float eps_) {
return copysign(1.0f, x_)*fmaxf(fabsf(x_), fminf(x_*x_/(2.0f*eps_)+0.5f*eps_, eps_));
}
/**
* Reads a block of data with ghost cells
*/
template<int vars, int sm_width, int sm_height, int block_width, int block_height>
inline __device__ void readBlock(float* ptr_[vars], int pitch_[vars],
float shmem[vars][sm_height][sm_width],
const int max_x_, const int max_y_) {
//Index of block within domain
const int bx = blockDim.x * blockIdx.x;
const int by = blockDim.y * blockIdx.y;
float* rows[3];
//Read into shared memory
for (int j=threadIdx.y; j<sm_height; j+=block_height) {
const int l = clamp(by + j, 0, max_y_-1); // Clamp out of bounds
//Compute the pointer to current row in the arrays
for (int m=0; m<vars; ++m) {
rows[m] = (float*) ((char*) ptr_[m] + pitch_[m]*l);
}
for (int i=threadIdx.x; i<sm_width; i+=block_width) {
const int k = clamp(bx + i, 0, max_x_-1); // Clamp out of bounds
for (int m=0; m<vars; ++m) {
shmem[m][j][i] = rows[m][k];
}
}
}
}
/**
* Writes a block of data to global memory for the shallow water equations.
*/
template<int sm_width, int sm_height, int offset_x=0, int offset_y=0>
inline __device__ void writeBlock(float* ptr_, int pitch_,
float shmem[sm_height][sm_width],
const int width, const int height) {
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + offset_x;
const int tj = blockDim.y*blockIdx.y + threadIdx.y + offset_y;
//Only write internal cells
if (ti < width+offset_x && tj < height+offset_y) {
//Index of thread within block
const int tx = threadIdx.x + offset_x;
const int ty = threadIdx.y + offset_y;
float* const row = (float*) ((char*) ptr_ + pitch_*tj);
row[ti] = shmem[ty][tx];
}
}
/**
* Writes a block of data to global memory for the shallow water equations.
*/
__device__ void writeBlock2(float* h_ptr_, int h_pitch_,
float* hu_ptr_, int hu_pitch_,
float* hv_ptr_, int hv_pitch_,
float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
const int nx_, const int ny_) {
writeBlock<BLOCK_WIDTH+4, BLOCK_HEIGHT+4, 2, 2>( h_ptr_, h_pitch_, Q[0], nx_, ny_);
writeBlock<BLOCK_WIDTH+4, BLOCK_HEIGHT+4, 2, 2>(hu_ptr_, hu_pitch_, Q[1], nx_, ny_);
writeBlock<BLOCK_WIDTH+4, BLOCK_HEIGHT+4, 2, 2>(hv_ptr_, hv_pitch_, Q[2], nx_, ny_);
}
/**
* No flow boundary conditions for the shallow water equations
* with one ghost cell in each direction
*/
__device__ void noFlowBoundary1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2], const int nx_, const int ny_) {
//Global index
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1;
//Block-local indices
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int i = tx + 1; //Skip local ghost cells, i.e., +1
const int j = ty + 1;
//Fix boundary conditions
if (ti == 1) {
Q[0][j][i-1] = Q[0][j][i];
Q[1][j][i-1] = -Q[1][j][i];
Q[2][j][i-1] = Q[2][j][i];
}
if (ti == nx_) {
Q[0][j][i+1] = Q[0][j][i];
Q[1][j][i+1] = -Q[1][j][i];
Q[2][j][i+1] = Q[2][j][i];
}
if (tj == 1) {
Q[0][j-1][i] = Q[0][j][i];
Q[1][j-1][i] = Q[1][j][i];
Q[2][j-1][i] = -Q[2][j][i];
}
if (tj == ny_) {
Q[0][j+1][i] = Q[0][j][i];
Q[1][j+1][i] = Q[1][j][i];
Q[2][j+1][i] = -Q[2][j][i];
}
}
/**
* No flow boundary conditions for the shallow water equations
* with two ghost cells in each direction
*/
__device__ void noFlowBoundary2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], const int nx_, const int ny_) {
//Global index
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2;
//Block-local indices
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int i = tx + 2; //Skip local ghost cells, i.e., +2
const int j = ty + 2;
if (ti == 2) {
Q[0][j][i-1] = Q[0][j][i];
Q[1][j][i-1] = -Q[1][j][i];
Q[2][j][i-1] = Q[2][j][i];
Q[0][j][i-2] = Q[0][j][i+1];
Q[1][j][i-2] = -Q[1][j][i+1];
Q[2][j][i-2] = Q[2][j][i+1];
}
if (ti == nx_+1) {
Q[0][j][i+1] = Q[0][j][i];
Q[1][j][i+1] = -Q[1][j][i];
Q[2][j][i+1] = Q[2][j][i];
Q[0][j][i+2] = Q[0][j][i-1];
Q[1][j][i+2] = -Q[1][j][i-1];
Q[2][j][i+2] = Q[2][j][i-1];
}
if (tj == 2) {
Q[0][j-1][i] = Q[0][j][i];
Q[1][j-1][i] = Q[1][j][i];
Q[2][j-1][i] = -Q[2][j][i];
Q[0][j-2][i] = Q[0][j+1][i];
Q[1][j-2][i] = Q[1][j+1][i];
Q[2][j-2][i] = -Q[2][j+1][i];
}
if (tj == ny_+1) {
Q[0][j+1][i] = Q[0][j][i];
Q[1][j+1][i] = Q[1][j][i];
Q[2][j+1][i] = -Q[2][j][i];
Q[0][j+2][i] = Q[0][j-1][i];
Q[1][j+2][i] = Q[1][j-1][i];
Q[2][j+2][i] = -Q[2][j-1][i];
}
}
/**
* Evolves the solution in time along the x axis (dimensional splitting)
*/
__device__ void evolveF1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
const int nx_, const int ny_,
const float dx_, const float dt_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1;
if (ti > 0 && ti < nx_+1 && tj > 0 && tj < ny_+1) {
const int i = tx + 1; //Skip local ghost cells, i.e., +1
const int j = ty + 1;
Q[0][j][i] = Q[0][j][i] + (F[0][ty][tx] - F[0][ty][tx+1]) * dt_ / dx_;
Q[1][j][i] = Q[1][j][i] + (F[1][ty][tx] - F[1][ty][tx+1]) * dt_ / dx_;
Q[2][j][i] = Q[2][j][i] + (F[2][ty][tx] - F[2][ty][tx+1]) * dt_ / dx_;
}
}
/**
* Evolves the solution in time along the x axis (dimensional splitting)
*/
__device__ void evolveF2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
const int nx_, const int ny_,
const float dx_, const float dt_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2;
if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) {
const int i = tx + 2; //Skip local ghost cells, i.e., +1
const int j = ty + 2;
Q[0][j][i] = Q[0][j][i] + (F[0][ty][tx] - F[0][ty][tx+1]) * dt_ / dx_;
Q[1][j][i] = Q[1][j][i] + (F[1][ty][tx] - F[1][ty][tx+1]) * dt_ / dx_;
Q[2][j][i] = Q[2][j][i] + (F[2][ty][tx] - F[2][ty][tx+1]) * dt_ / dx_;
}
}
/**
* Evolves the solution in time along the y axis (dimensional splitting)
*/
__device__ void evolveG1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
const int nx_, const int ny_,
const float dy_, const float dt_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1;
if (ti > 0 && ti < nx_+1 && tj > 0 && tj < ny_+1) {
const int i = tx + 1; //Skip local ghost cells, i.e., +1
const int j = ty + 1;
Q[0][j][i] = Q[0][j][i] + (G[0][ty][tx] - G[0][ty+1][tx]) * dt_ / dy_;
Q[1][j][i] = Q[1][j][i] + (G[1][ty][tx] - G[1][ty+1][tx]) * dt_ / dy_;
Q[2][j][i] = Q[2][j][i] + (G[2][ty][tx] - G[2][ty+1][tx]) * dt_ / dy_;
}
}
/**
* Evolves the solution in time along the y axis (dimensional splitting)
*/
__device__ void evolveG2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
const int nx_, const int ny_,
const float dy_, const float dt_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2
const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2;
if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) {
const int i = tx + 2; //Skip local ghost cells, i.e., +2
const int j = ty + 2;
Q[0][j][i] = Q[0][j][i] + (G[0][ty][tx] - G[0][ty+1][tx]) * dt_ / dy_;
Q[1][j][i] = Q[1][j][i] + (G[1][ty][tx] - G[1][ty+1][tx]) * dt_ / dy_;
Q[2][j][i] = Q[2][j][i] + (G[2][ty][tx] - G[2][ty+1][tx]) * dt_ / dy_;
}
}

129
GPUSimulators/cuda/limiters.h Normal file
View File

@@ -0,0 +1,129 @@
/*
This file implements different flux and slope limiters
Copyright (C) 2016, 2017, 2018 SINTEF ICT
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/>.
*/
#pragma once
/**
* Reconstructs a slope using the generalized minmod limiter based on three
* consecutive values
*/
__device__ __inline__ float minmodSlope(float left, float center, float right, float theta) {
const float backward = (center - left) * theta;
const float central = (right - left) * 0.5f;
const float forward = (right - center) * theta;
return 0.25f
*copysign(1.0f, backward)
*(copysign(1.0f, backward) + copysign(1.0f, central))
*(copysign(1.0f, central) + copysign(1.0f, forward))
*min( min(fabs(backward), fabs(central)), fabs(forward) );
}
/**
* Reconstructs a minmod slope for a whole block along the abscissa
*/
__device__ void minmodSlopeX(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
const float theta_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
//Reconstruct slopes along x axis
{
const int j = ty;
const int l = j + 2; //Skip ghost cells
for (int i=tx; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
const int k = i + 1;
for (int p=0; p<3; ++p) {
Qx[p][j][i] = minmodSlope(Q[p][l][k-1], Q[p][l][k], Q[p][l][k+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_) {
//Index of thread within block
const int tx = threadIdx.x;
const int ty = threadIdx.y;
for (int j=ty; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
const int l = j + 1;
{
const int i = tx;
const int k = i + 2; //Skip ghost cells
for (int p=0; p<3; ++p) {
Qy[p][j][i] = minmodSlope(Q[p][l-1][k], Q[p][l][k], Q[p][l+1][k], theta_);
}
}
}
}
__device__ float monotonized_central(float r_) {
return fmaxf(0.0f, fminf(2.0f, fminf(2.0f*r_, 0.5f*(1.0f+r_))));
}
__device__ float osher(float r_, float beta_) {
return fmaxf(0.0f, fminf(beta_, r_));
}
__device__ float sweby(float r_, float beta_) {
return fmaxf(0.0f, fmaxf(fminf(r_, beta_), fminf(beta_*r_, 1.0f)));
}
__device__ float minmod(float r_) {
return fmaxf(0.0f, fminf(1.0f, r_));
}
__device__ float generalized_minmod(float r_, float theta_) {
return fmaxf(0.0f, fminf(theta_*r_, fminf( (1.0f + r_) / 2.0f, theta_)));
}
__device__ float superbee(float r_) {
return fmaxf(0.0f, fmaxf(fminf(2.0f*r_, 1.0f), fminf(r_, 2.0f)));
}
__device__ float vanAlbada1(float r_) {
return (r_*r_ + r_) / (r_*r_ + 1.0f);
}
__device__ float vanAlbada2(float r_) {
return 2.0f*r_ / (r_*r_* + 1.0f);
}
__device__ float vanLeer(float r_) {
return (r_ + fabsf(r_)) / (1.0f + fabsf(r_));
}