smyalygames/FiniteVolumeGPU
1
0
Fork 0
mirror of https://github.com/smyalygames/FiniteVolumeGPU.git synced 2025-05-18 06:24:13 +02:00
Code Issues Actions Packages Projects Releases Wiki Activity

Updated compilation of kernels etc

Browse Source
This commit is contained in:
André R. Brodtkorb 2018-07-25 15:06:56 +02:00
parent 7cecd23e59
commit c6758b477b
7 changed files with 511 additions and 513 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

729
ConvergenceSmooth1D.ipynb
View File

File diff suppressed because one or more lines are too long

103
SWESimulators/Common.py
View File

@ -1,40 +1,99 @@
import os
import numpy as np
import time
import pycuda.compiler as cuda_compiler
import pycuda.gpuarray
import pycuda.driver as cuda
"""
Static function which reads a text file and creates an OpenCL kernel from that
Class which keeps track of the CUDA context and some helper functions
"""
def get_kernel(kernel_filename, block_width, block_height):
import datetime
#Create define string
define_string = "#define block_width " + str(block_width) + "\n"
define_string += "#define block_height " + str(block_height) + "\n\n"
#Read the proper program
module_path = os.path.dirname(os.path.realpath(__file__))
fullpath = os.path.join(module_path, kernel_filename)
#with open(fullpath, "r") as kernel_file:
# kernel_string = define_string + kernel_file.read()
# kernel = cuda_compiler.SourceModule(kernel_string, include_dirs=[module_path], no_extern_c=True)
kernel_string = define_string + '#include "' + fullpath + '"'
kernel = cuda_compiler.SourceModule(kernel_string, include_dirs=[module_path])
class CudaContext(object):
def __init__(self, verbose=True, blocking=False):
self.verbose = verbose
self.blocking = blocking
self.kernels = {}
return kernel
cuda.init(flags=0)
if (self.verbose):
print("CUDA version " + str(cuda.get_version()))
print("Driver version " + str(cuda.get_driver_version()))
self.cuda_device = cuda.Device(0)
if (self.verbose):
print("Using " + self.cuda_device.name())
print(" => compute capability: " + str(self.cuda_device.compute_capability()))
print(" => memory: " + str(self.cuda_device.total_memory() / (1024*1024)) + " MB")
if (self.blocking):
self.cuda_context = self.cuda_device.make_context(flags=cuda.ctx_flags.SCHED_BLOCKING_SYNC)
if (self.verbose):
print("=== WARNING ===")
print("Using blocking context")
print("=== WARNING ===")
else:
self.cuda_context = self.cuda_device.make_context(flags=cuda.ctx_flags.SCHED_AUTO)
def __del__(self, *args):
if self.verbose:
print("Cleaning up CUDA context")
self.cuda_context.detach()
cuda.Context.pop()
"""
Reads a text file and creates an OpenCL kernel from that
"""
def get_kernel(self, kernel_filename, block_width, block_height):
# Generate a kernel ID for our cache
module_path = os.path.dirname(os.path.realpath(__file__))
fullpath = os.path.join(module_path, kernel_filename)
kernel_date = os.path.getmtime(fullpath)
with open(fullpath, "r") as kernel_file:
kernel_hash = hash(kernel_file.read())
kernel_id = kernel_filename + ":" + str(kernel_hash) + ":" + str(kernel_date)
# Simple caching to keep keep from recompiling kernels
if (kernel_id not in self.kernels.keys()):
#Create define string
define_string = "#define block_width " + str(block_width) + "\n"
define_string += "#define block_height " + str(block_height) + "\n\n"
kernel_string = define_string + '#include "' + fullpath + '"'
self.kernels[kernel_id] = cuda_compiler.SourceModule(kernel_string, include_dirs=[module_path])
return self.kernels[kernel_id]
"""
Clears the kernel cache (useful for debugging & development)
"""
def clear_kernel_cache(self):
self.kernels = {}
class Timer(object):
def __init__(self, tag, verbose=True):
self.verbose = verbose
self.tag = tag
def __enter__(self):
self.start = time.time()
return self
def __exit__(self, *args):
self.end = time.time()
self.secs = self.end - self.start
self.msecs = self.secs * 1000 # millisecs
if self.verbose:
print("=> " + self.tag + ' %f ms' % self.msecs)
@ -53,9 +112,9 @@ class CUDAArray2D:
self.ny_halo = ny + 2*halo_y
#Make sure data is in proper format
assert(np.issubdtype(data.dtype, np.float32), "Wrong datatype: %s" % str(data.dtype))
assert(not np.isfortran(data), "Wrong datatype (Fortran, expected C)")
assert(data.shape == (self.ny_halo, self.nx_halo), "Wrong data shape: %s" % str(data.shape))
assert np.issubdtype(data.dtype, np.float32), "Wrong datatype: %s" % str(data.dtype)
assert not np.isfortran(data), "Wrong datatype (Fortran, expected C)"
assert data.shape == (self.ny_halo, self.nx_halo), "Wrong data shape: %s" % str(data.shape)
#Upload data to the device
self.data = pycuda.gpuarray.to_gpu_async(data, stream=stream)

3
SWESimulators/FORCE.py
View File

@ -57,6 +57,7 @@ class FORCE:
g: Gravitational accelleration (9.81 m/s^2)
"""
def __init__(self, \
context, \
h0, hu0, hv0, \
nx, ny, \
dx, dy, dt, \
@ -66,7 +67,7 @@ class FORCE:
self.stream = cuda.Stream()
#Get kernels
self.force_module = Common.get_kernel("FORCE_kernel.cu", block_width, block_height)
self.force_module = context.get_kernel("FORCE_kernel.cu", block_width, block_height)
self.force_kernel = self.force_module.get_function("FORCEKernel")
self.force_kernel.prepare("iiffffPiPiPiPiPiPi")

3
SWESimulators/HLL.py
View File

@ -52,6 +52,7 @@ class HLL:
g: Gravitational accelleration (9.81 m/s^2)
"""
def __init__(self, \
context, \
h0, hu0, hv0, \
nx, ny, \
dx, dy, dt, \
@ -61,7 +62,7 @@ class HLL:
self.stream = cuda.Stream()
#Get kernels
self.hll_module = Common.get_kernel("HLL_kernel.cu", block_width, block_height)
self.hll_module = context.get_kernel("HLL_kernel.cu", block_width, block_height)
self.hll_kernel = self.hll_module.get_function("HLLKernel")
self.hll_kernel.prepare("iiffffPiPiPiPiPiPi")

67
SWESimulators/HLL2.py
View File

@ -21,7 +21,11 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#Import packages we need
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
@ -50,25 +54,28 @@ class HLL2:
g: Gravitational accelleration (9.81 m/s^2)
"""
def __init__(self, \
cl_ctx, \
context, \
h0, hu0, hv0, \
nx, ny, \
dx, dy, dt, \
g, \
theta=1.8, \
block_width=16, block_height=16):
self.cl_ctx = cl_ctx
#Create an OpenCL command queue
self.cl_queue = cl.CommandQueue(self.cl_ctx)
#Create a CUDA stream
self.stream = cuda.Stream()
#Get kernels
self.swe_kernel = Common.get_kernel(self.cl_ctx, "HLL2_kernel.opencl", block_width, block_height)
self.hll2_module = context.get_kernel("HLL2_kernel.cu", block_width, block_height)
self.hll2_kernel = self.hll2_module.get_function("HLL2Kernel")
self.hll2_kernel.prepare("iifffffiPiPiPiPiPiPi")
#Create data by uploading to device
ghost_cells_x = 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
#Notice that we need to specify them in the correct dataformat for the
@ -85,15 +92,15 @@ class HLL2:
self.t = np.float32(0.0)
#Compute kernel launch parameters
self.local_size = (block_width, block_height)
self.local_size = (block_width, block_height, 1)
self.global_size = ( \
int(np.ceil(self.nx / float(self.local_size[0])) * self.local_size[0]), \
int(np.ceil(self.ny / float(self.local_size[1])) * self.local_size[1]) \
)
int(np.ceil(self.nx / float(self.local_size[0]))), \
int(np.ceil(self.ny / float(self.local_size[1]))) \
)
def __str__(self):
return "Harten-Lax-van Leer contact discontinuity"
return "Harten-Lax-van Leer (2nd order)"
"""
@ -111,34 +118,34 @@ class HLL2:
break
#Along X, then Y
self.swe_kernel.swe_2D(self.cl_queue, self.global_size, self.local_size, \
self.hll2_kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
self.nx, self.ny, \
self.dx, self.dy, local_dt, \
self.g, \
self.theta, \
np.int32(0), \
self.cl_data.h0.data, self.cl_data.h0.pitch, \
self.cl_data.hu0.data, self.cl_data.hu0.pitch, \
self.cl_data.hv0.data, self.cl_data.hv0.pitch, \
self.cl_data.h1.data, self.cl_data.h1.pitch, \
self.cl_data.hu1.data, self.cl_data.hu1.pitch, \
self.cl_data.hv1.data, self.cl_data.hv1.pitch)
self.cl_data.swap()
self.data.h0.data.gpudata, self.data.h0.pitch, \
self.data.hu0.data.gpudata, self.data.hu0.pitch, \
self.data.hv0.data.gpudata, self.data.hv0.pitch, \
self.data.h1.data.gpudata, self.data.h1.pitch, \
self.data.hu1.data.gpudata, self.data.hu1.pitch, \
self.data.hv1.data.gpudata, self.data.hv1.pitch)
self.data.swap()
#Along Y, then X
self.swe_kernel.swe_2D(self.cl_queue, self.global_size, self.local_size, \
self.hll2_kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
self.nx, self.ny, \
self.dx, self.dy, local_dt, \
self.g, \
self.theta, \
np.int32(1), \
self.cl_data.h0.data, self.cl_data.h0.pitch, \
self.cl_data.hu0.data, self.cl_data.hu0.pitch, \
self.cl_data.hv0.data, self.cl_data.hv0.pitch, \
self.cl_data.h1.data, self.cl_data.h1.pitch, \
self.cl_data.hu1.data, self.cl_data.hu1.pitch, \
self.cl_data.hv1.data, self.cl_data.hv1.pitch)
self.cl_data.swap()
self.data.h0.data.gpudata, self.data.h0.pitch, \
self.data.hu0.data.gpudata, self.data.hu0.pitch, \
self.data.hv0.data.gpudata, self.data.hv0.pitch, \
self.data.h1.data.gpudata, self.data.h1.pitch, \
self.data.hu1.data.gpudata, self.data.hu1.pitch, \
self.data.hv1.data.gpudata, self.data.hv1.pitch)
self.data.swap()
self.t += local_dt
@ -148,5 +155,5 @@ class HLL2:
def download(self):
return self.cl_data.download(self.cl_queue)
return self.data.download(self.stream)

116
SWESimulators/HLL2_kernel.opencl → SWESimulators/HLL2_kernel.cu
View File

@ -18,7 +18,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "common.opencl"
#include "common.cu"
@ -29,9 +29,10 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
/**
* Computes the flux along the x axis for all faces
*/
void computeFluxF(__local float Q[3][block_height+4][block_width+4],
__local float Qx[3][block_height+2][block_width+2],
__local float F[3][block_height+1][block_width+1],
__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 = get_local_id(0);
@ -43,19 +44,19 @@ void computeFluxF(__local float Q[3][block_height+4][block_width+4],
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 = (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 = (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][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 = (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 = (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][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_));
@ -79,9 +80,10 @@ void computeFluxF(__local float Q[3][block_height+4][block_width+4],
/**
* Computes the flux along the x axis for all faces
*/
void computeFluxG(__local float Q[3][block_height+4][block_width+4],
__local float Qy[3][block_height+2][block_width+2],
__local float G[3][block_height+1][block_width+1],
__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 = get_local_id(0);
@ -94,19 +96,19 @@ void computeFluxG(__local float Q[3][block_height+4][block_width+4],
// 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 = (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 = (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][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 = (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 = (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][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_));
@ -131,7 +133,7 @@ void computeFluxG(__local float Q[3][block_height+4][block_width+4],
__kernel void swe_2D(
__global__ void HLL2Kernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
@ -141,19 +143,19 @@ __kernel void swe_2D(
int step_,
//Input h^n
__global float* h0_ptr_, int h0_pitch_,
__global float* hu0_ptr_, int hu0_pitch_,
__global float* hv0_ptr_, int hv0_pitch_,
float* h0_ptr_, int h0_pitch_,
float* hu0_ptr_, int hu0_pitch_,
float* hv0_ptr_, int hv0_pitch_,
//Output h^{n+1}
__global float* h1_ptr_, int h1_pitch_,
__global float* hu1_ptr_, int hu1_pitch_,
__global float* hv1_ptr_, int hv1_pitch_) {
float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) {
//Shared memory variables
__local float Q[3][block_height+4][block_width+4];
__local float Qx[3][block_height+2][block_width+2];
__local float F[3][block_height+1][block_width+1];
__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];
@ -163,55 +165,55 @@ __kernel void swe_2D(
hu0_ptr_, hu0_pitch_,
hv0_ptr_, hv0_pitch_,
Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Step 0 => evolve x first, then y
if (step_ == 0) {
//Compute fluxes along the x axis and evolve
minmodSlopeX(Q, Qx, theta_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Compute fluxes along the y axis and evolve
minmodSlopeY(Q, Qx, theta_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
}
//Step 1 => evolve y first, then x
else {
//Compute fluxes along the y axis and evolve
minmodSlopeY(Q, Qx, theta_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
evolveG2(Q, F, nx_, ny_, dy_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Set boundary conditions
noFlowBoundary2(Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Compute fluxes along the x axis and evolve
minmodSlopeX(Q, Qx, theta_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
evolveF2(Q, F, nx_, ny_, dx_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
}

3
SWESimulators/LxF.py
View File

@ -53,6 +53,7 @@ class LxF:
g: Gravitational accelleration (9.81 m/s^2)
"""
def __init__(self, \
context, \
h0, hu0, hv0, \
nx, ny, \
dx, dy, dt, \
@ -62,7 +63,7 @@ class LxF:
self.stream = None #cuda.Stream()
#Get kernels
self.lxf_module = Common.get_kernel("LxF_kernel.cu", block_width, block_height)
self.lxf_module = context.get_kernel("LxF_kernel.cu", block_width, block_height)
self.lxf_kernel = self.lxf_module.get_function("LxFKernel")
self.lxf_kernel.prepare("iiffffPiPiPiPiPiPi")