smyalygames/FiniteVolumeGPU
1
0
Fork 0
mirror of https://github.com/smyalygames/FiniteVolumeGPU.git synced 2025-11-26 08:09:49 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity

Fixed KP07

Browse Source
This commit is contained in:
André R. Brodtkorb 2018-07-25 15:52:31 +02:00
parent c6758b477b
commit cbb1bdb839
4 changed files with 277 additions and 273 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

338
ConvergenceSmooth1D.ipynb
View File

File diff suppressed because one or more lines are too long

98
SWESimulators/KP07.py
View File

@ -26,7 +26,11 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#Import packages we need #Import packages we need
import numpy as np import numpy as np
import pyopencl as cl #OpenCL in Python
import pycuda.compiler as cuda_compiler
import pycuda.gpuarray
import pycuda.driver as cuda
from SWESimulators import Common from SWESimulators import Common
@ -63,26 +67,25 @@ class KP07:
wind_v0: Translation speed along y for moving cyclone (-0.5*u0) wind_v0: Translation speed along y for moving cyclone (-0.5*u0)
""" """
def __init__(self, \ def __init__(self, \
cl_ctx, \ context, \
h0, hu0, hv0, \ h0, hu0, hv0, \
nx, ny, \ nx, ny, \
dx, dy, dt, \ dx, dy, dt, \
g, f=0.0, r=0.0, \ g, theta=1.3, \
theta=1.3, use_rk2=True, r=0.0, use_rk2=True,
wind_stress=Common.WindStressParams(), \
block_width=16, block_height=16): block_width=16, block_height=16):
self.cl_ctx = cl_ctx #Create a CUDA stream
self.stream = cuda.Stream()
#Create an OpenCL command queue
self.cl_queue = cl.CommandQueue(self.cl_ctx)
#Get kernels #Get kernels
self.kp07_kernel = Common.get_kernel(self.cl_ctx, "KP07_kernel.opencl", block_width, block_height) self.kp07_module = context.get_kernel("KP07_kernel.cu", block_width, block_height)
self.kp07_kernel = self.kp07_module.get_function("KP07Kernel")
self.kp07_kernel.prepare("iiffffffiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
ghost_cells_x = 2 ghost_cells_x = 2
ghost_cells_y = 2 ghost_cells_y = 2
self.cl_data = Common.SWEDataArkawaA(self.cl_ctx, nx, ny, ghost_cells_x, ghost_cells_y, h0, hu0, hv0) self.data = Common.SWEDataArakawaA(self.stream, nx, ny, ghost_cells_x, ghost_cells_y, h0, hu0, hv0)
#Save input parameters #Save input parameters
#Notice that we need to specify them in the correct dataformat for the #Notice that we need to specify them in the correct dataformat for the
@ -93,26 +96,24 @@ class KP07:
self.dy = np.float32(dy) self.dy = np.float32(dy)
self.dt = np.float32(dt) self.dt = np.float32(dt)
self.g = np.float32(g) self.g = np.float32(g)
self.f = np.float32(f)
self.r = np.float32(r) self.r = np.float32(r)
self.theta = np.float32(theta) self.theta = np.float32(theta)
self.use_rk2 = use_rk2 self.use_rk2 = use_rk2
self.wind_stress = wind_stress
#Initialize time #Initialize time
self.t = np.float32(0.0) self.t = np.float32(0.0)
#Compute kernel launch parameters #Compute kernel launch parameters
self.local_size = (block_width, block_height) self.local_size = (block_width, block_height, 1)
self.global_size = ( \ self.global_size = ( \
int(np.ceil(self.nx / float(self.local_size[0])) * self.local_size[0]), \ int(np.ceil(self.nx / float(self.local_size[0]))), \
int(np.ceil(self.ny / float(self.local_size[1])) * self.local_size[1]) \ int(np.ceil(self.ny / float(self.local_size[1]))) \
) )
def __str__(self): def __str__(self):
return "Kurganov-Petrova" return "Kurganov-Petrova 2007"
""" """
Function which steps n timesteps Function which steps n timesteps
@ -127,64 +128,47 @@ class KP07:
break break
if (self.use_rk2): if (self.use_rk2):
self.kp07_kernel.swe_2D(self.cl_queue, self.global_size, self.local_size, \ self.kp07_kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
self.nx, self.ny, \ self.nx, self.ny, \
self.dx, self.dy, local_dt, \ self.dx, self.dy, local_dt, \
self.g, \ self.g, \
self.theta, \ self.theta, \
self.f, \
self.r, \ self.r, \
np.int32(0), \ np.int32(0), \
self.cl_data.h0.data, self.cl_data.h0.pitch, \ self.data.h0.data.gpudata, self.data.h0.pitch, \
self.cl_data.hu0.data, self.cl_data.hu0.pitch, \ self.data.hu0.data.gpudata, self.data.hu0.pitch, \
self.cl_data.hv0.data, self.cl_data.hv0.pitch, \ self.data.hv0.data.gpudata, self.data.hv0.pitch, \
self.cl_data.h1.data, self.cl_data.h1.pitch, \ self.data.h1.data.gpudata, self.data.h1.pitch, \
self.cl_data.hu1.data, self.cl_data.hu1.pitch, \ self.data.hu1.data.gpudata, self.data.hu1.pitch, \
self.cl_data.hv1.data, self.cl_data.hv1.pitch, \ self.data.hv1.data.gpudata, self.data.hv1.pitch)
self.wind_stress.type, \
self.wind_stress.tau0, self.wind_stress.rho, self.wind_stress.alpha, self.wind_stress.xm, self.wind_stress.Rc, \ self.kp07_kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
self.wind_stress.x0, self.wind_stress.y0, \
self.wind_stress.u0, self.wind_stress.v0, \
self.t)
self.kp07_kernel.swe_2D(self.cl_queue, self.global_size, self.local_size, \
self.nx, self.ny, \ self.nx, self.ny, \
self.dx, self.dy, local_dt, \ self.dx, self.dy, local_dt, \
self.g, \ self.g, \
self.theta, \ self.theta, \
self.f, \
self.r, \ self.r, \
np.int32(1), \ np.int32(1), \
self.cl_data.h1.data, self.cl_data.h1.pitch, \ self.data.h1.data.gpudata, self.data.h1.pitch, \
self.cl_data.hu1.data, self.cl_data.hu1.pitch, \ self.data.hu1.data.gpudata, self.data.hu1.pitch, \
self.cl_data.hv1.data, self.cl_data.hv1.pitch, \ self.data.hv1.data.gpudata, self.data.hv1.pitch, \
self.cl_data.h0.data, self.cl_data.h0.pitch, \ self.data.h0.data.gpudata, self.data.h0.pitch, \
self.cl_data.hu0.data, self.cl_data.hu0.pitch, \ self.data.hu0.data.gpudata, self.data.hu0.pitch, \
self.cl_data.hv0.data, self.cl_data.hv0.pitch, \ self.data.hv0.data.gpudata, self.data.hv0.pitch)
self.wind_stress.type, \
self.wind_stress.tau0, self.wind_stress.rho, self.wind_stress.alpha, self.wind_stress.xm, self.wind_stress.Rc, \
self.wind_stress.x0, self.wind_stress.y0, \
self.wind_stress.u0, self.wind_stress.v0, \
self.t)
else: else:
self.kp07_kernel.swe_2D(self.cl_queue, self.global_size, self.local_size, \ self.kp07_kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
self.nx, self.ny, \ self.nx, self.ny, \
self.dx, self.dy, local_dt, \ self.dx, self.dy, local_dt, \
self.g, \ self.g, \
self.theta, \ self.theta, \
self.f, \
self.r, \ self.r, \
np.int32(0), \ np.int32(0), \
self.cl_data.h0.data, self.cl_data.h0.pitch, \ self.data.h0.data.gpudata, self.data.h0.pitch, \
self.cl_data.hu0.data, self.cl_data.hu0.pitch, \ self.data.hu0.data.gpudata, self.data.hu0.pitch, \
self.cl_data.hv0.data, self.cl_data.hv0.pitch, \ self.data.hv0.data.gpudata, self.data.hv0.pitch, \
self.cl_data.h1.data, self.cl_data.h1.pitch, \ self.data.h1.data.gpudata, self.data.h1.pitch, \
self.cl_data.hu1.data, self.cl_data.hu1.pitch, \ self.data.hu1.data.gpudata, self.data.hu1.pitch, \
self.cl_data.hv1.data, self.cl_data.hv1.pitch, \ self.data.hv1.data.gpudata, self.data.hv1.pitch)
self.wind_stress.type, \
self.wind_stress.tau0, self.wind_stress.rho, self.wind_stress.alpha, self.wind_stress.xm, self.wind_stress.Rc, \
self.wind_stress.x0, self.wind_stress.y0, \
self.wind_stress.u0, self.wind_stress.v0, \
self.t)
self.cl_data.swap() self.cl_data.swap()
self.t += local_dt self.t += local_dt
@ -196,5 +180,5 @@ class KP07:
def download(self): def download(self):
return self.cl_data.download(self.cl_queue) return self.data.download(self.stream)

96
SWESimulators/KP07_kernel.opencl → SWESimulators/KP07_kernel.cu
View File

@ -24,13 +24,13 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "common.opencl" #include "common.cu"
__device__
void computeFluxF(__local float Q[3][block_height+4][block_width+4], void computeFluxF(float Q[3][block_height+4][block_width+4],
__local float Qx[3][block_height+2][block_width+2], float Qx[3][block_height+2][block_width+2],
__local float F[3][block_height+1][block_width+1], float F[3][block_height+1][block_width+1],
const float g_) { const float g_) {
//Index of thread within block //Index of thread within block
const int tx = get_local_id(0); const int tx = get_local_id(0);
@ -41,10 +41,10 @@ void computeFluxF(__local float Q[3][block_height+4][block_width+4],
for (int i=tx; i<block_width+1; i+=get_local_size(0)) { for (int i=tx; i<block_width+1; i+=get_local_size(0)) {
const int k = i + 1; const int k = i + 1;
// Q at interface from the right and left // Q at interface from the right and left
const float3 Qp = (float3)(Q[0][l][k+1] - 0.5f*Qx[0][j][i+1], 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[1][l][k+1] - 0.5f*Qx[1][j][i+1],
Q[2][l][k+1] - 0.5f*Qx[2][j][i+1]); Q[2][l][k+1] - 0.5f*Qx[2][j][i+1]);
const float3 Qm = (float3)(Q[0][l][k ] + 0.5f*Qx[0][j][i ], 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[1][l][k ] + 0.5f*Qx[1][j][i ],
Q[2][l][k ] + 0.5f*Qx[2][j][i ]); Q[2][l][k ] + 0.5f*Qx[2][j][i ]);
@ -57,9 +57,10 @@ void computeFluxF(__local float Q[3][block_height+4][block_width+4],
} }
} }
void computeFluxG(__local float Q[3][block_height+4][block_width+4], __device__
__local float Qy[3][block_height+2][block_width+2], void computeFluxG(float Q[3][block_height+4][block_width+4],
__local float G[3][block_height+1][block_width+1], float Qy[3][block_height+2][block_width+2],
float G[3][block_height+1][block_width+1],
const float g_) { const float g_) {
//Index of thread within block //Index of thread within block
const int tx = get_local_id(0); const int tx = get_local_id(0);
@ -71,10 +72,10 @@ void computeFluxG(__local float Q[3][block_height+4][block_width+4],
const int k = i + 2; //Skip ghost cells const int k = i + 2; //Skip ghost cells
// Q at interface from the right and left // Q at interface from the right and left
// Note that we swap hu and hv // Note that we swap hu and hv
const float3 Qp = (float3)(Q[0][l+1][k] - 0.5f*Qy[0][j+1][i], 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[2][l+1][k] - 0.5f*Qy[2][j+1][i],
Q[1][l+1][k] - 0.5f*Qy[1][j+1][i]); Q[1][l+1][k] - 0.5f*Qy[1][j+1][i]);
const float3 Qm = (float3)(Q[0][l ][k] + 0.5f*Qy[0][j ][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[2][l ][k] + 0.5f*Qy[2][j ][i],
Q[1][l ][k] + 0.5f*Qy[1][j ][i]); Q[1][l ][k] + 0.5f*Qy[1][j ][i]);
@ -94,56 +95,44 @@ void computeFluxG(__local float Q[3][block_height+4][block_width+4],
/** /**
* This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme * This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
*/ */
__kernel void swe_2D( __global__ void KP07Kernel(
int nx_, int ny_, int nx_, int ny_,
float dx_, float dy_, float dt_, float dx_, float dy_, float dt_,
float g_, float g_,
float theta_, float theta_,
float f_, //< Coriolis coefficient
float r_, //< Bottom friction coefficient float r_, //< Bottom friction coefficient
int step_, int step_,
//Input h^n //Input h^n
__global float* h0_ptr_, int h0_pitch_, float* h0_ptr_, int h0_pitch_,
__global float* hu0_ptr_, int hu0_pitch_, float* hu0_ptr_, int hu0_pitch_,
__global float* hv0_ptr_, int hv0_pitch_, float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1} //Output h^{n+1}
__global float* h1_ptr_, int h1_pitch_, float* h1_ptr_, int h1_pitch_,
__global float* hu1_ptr_, int hu1_pitch_, float* hu1_ptr_, int hu1_pitch_,
__global float* hv1_ptr_, int hv1_pitch_, float* hv1_ptr_, int hv1_pitch_) {
//Wind stress parameters
int wind_stress_type_,
float tau0_, float rho_, float alpha_, float xm_, float Rc_,
float x0_, float y0_,
float u0_, float v0_,
float t_) {
//Index of thread within block //Index of thread within block
const int tx = get_local_id(0); const int tx = get_local_id(0);
const int ty = get_local_id(1); const int ty = get_local_id(1);
//Index of block within domain
const int bx = get_local_size(0) * get_group_id(0);
const int by = get_local_size(1) * get_group_id(1);
//Index of cell within domain //Index of cell within domain
const int ti = get_global_id(0) + 2; //Skip global ghost cells, i.e., +2 const int ti = get_global_id(0) + 2; //Skip global ghost cells, i.e., +2
const int tj = get_global_id(1) + 2; const int tj = get_global_id(1) + 2;
//Shared memory variables //Shared memory variables
__local float Q[3][block_height+4][block_width+4]; __shared__ float Q[3][block_height+4][block_width+4];
//The following slightly wastes memory, but enables us to reuse the //The following slightly wastes memory, but enables us to reuse the
//funcitons in common.opencl //funcitons in common.opencl
__local float Qx[3][block_height+2][block_width+2]; __shared__ float Qx[3][block_height+2][block_width+2];
__local float Qy[3][block_height+2][block_width+2]; __shared__ float Qy[3][block_height+2][block_width+2];
__local float F[3][block_height+1][block_width+1]; __shared__ float F[3][block_height+1][block_width+1];
__local float G[3][block_height+1][block_width+1]; __shared__ float G[3][block_height+1][block_width+1];
@ -152,24 +141,24 @@ __kernel void swe_2D(
hu0_ptr_, hu0_pitch_, hu0_ptr_, hu0_pitch_,
hv0_ptr_, hv0_pitch_, hv0_ptr_, hv0_pitch_,
Q, nx_, ny_); Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE); __syncthreads();
//Fix boundary conditions //Fix boundary conditions
noFlowBoundary2(Q, nx_, ny_); noFlowBoundary2(Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE); __syncthreads();
//Reconstruct slopes along x and axis //Reconstruct slopes along x and axis
minmodSlopeX(Q, Qx, theta_); minmodSlopeX(Q, Qx, theta_);
minmodSlopeY(Q, Qy, theta_); minmodSlopeY(Q, Qy, theta_);
barrier(CLK_LOCAL_MEM_FENCE); __syncthreads();
//Compute fluxes along the x and y axis //Compute fluxes along the x and y axis
computeFluxF(Q, Qx, F, g_); computeFluxF(Q, Qx, F, g_);
computeFluxG(Q, Qy, G, g_); computeFluxG(Q, Qy, G, g_);
barrier(CLK_LOCAL_MEM_FENCE); __syncthreads();
//Sum fluxes and advance in time for all internal cells //Sum fluxes and advance in time for all internal cells
@ -177,33 +166,16 @@ __kernel void swe_2D(
const int i = tx + 2; //Skip local ghost cells, i.e., +2 const int i = tx + 2; //Skip local ghost cells, i.e., +2
const int j = ty + 2; const int j = ty + 2;
const float X = windStressX(
wind_stress_type_,
dx_, dy_, dt_,
tau0_, rho_, alpha_, xm_, Rc_,
x0_, y0_,
u0_, v0_,
t_);
const float Y = windStressY(
wind_stress_type_,
dx_, dy_, dt_,
tau0_, rho_, alpha_, xm_, Rc_,
x0_, y0_,
u0_, v0_,
t_);
const float h1 = Q[0][j][i] + (F[0][ty][tx] - F[0][ty ][tx+1]) * dt_ / dx_ 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_; + (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_ 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_ + (G[1][ty][tx] - G[1][ty+1][tx ]) * dt_ / dy_;
+ dt_*X - dt_*f_*Q[2][j][i];
const float hv1 = Q[2][j][i] + (F[2][ty][tx] - F[2][ty ][tx+1]) * dt_ / dx_ 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_ + (G[2][ty][tx] - G[2][ty+1][tx ]) * dt_ / dy_;
+ dt_*Y + dt_*f_*Q[1][j][i];
__global float* const h_row = (__global float*) ((__global char*) h1_ptr_ + h1_pitch_*tj); float* const h_row = (float*) ((char*) h1_ptr_ + h1_pitch_*tj);
__global float* const hu_row = (__global float*) ((__global char*) hu1_ptr_ + hu1_pitch_*tj); float* const hu_row = (float*) ((char*) hu1_ptr_ + hu1_pitch_*tj);
__global float* const hv_row = (__global float*) ((__global char*) hv1_ptr_ + hv1_pitch_*tj); float* const hv_row = (float*) ((char*) hv1_ptr_ + hv1_pitch_*tj);
const float C = 2.0f*r_*dt_/Q[0][j][i]; const float C = 2.0f*r_*dt_/Q[0][j][i];

2
SWESimulators/LxF.py
View File

@ -60,7 +60,7 @@ class LxF:
g, \ g, \
block_width=16, block_height=16): block_width=16, block_height=16):
#Create a CUDA stream #Create a CUDA stream
self.stream = None #cuda.Stream() self.stream = cuda.Stream()
#Get kernels #Get kernels
self.lxf_module = context.get_kernel("LxF_kernel.cu", block_width, block_height) self.lxf_module = context.get_kernel("LxF_kernel.cu", block_width, block_height)