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

Ported FORCE to CUDA

Browse Source
This commit is contained in:
André R. Brodtkorb 2018-07-25 08:48:43 +02:00
parent bc086865de
commit fcc1d0db1c
3 changed files with 114 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
OpenCL to CUDA.ipynb
View File

File diff suppressed because one or more lines are too long

42
SWESimulators/FORCE.py
View File

@ -22,7 +22,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
@ -53,24 +57,26 @@ class FORCE:
g: Gravitational accelleration (9.81 m/s^2)
"""
def __init__(self, \
cl_ctx, \
h0, hu0, hv0, \
nx, ny, \
dx, dy, dt, \
g, \
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.kernel = Common.get_kernel(self.cl_ctx, "FORCE_kernel.opencl", block_width, block_height)
self.force_module = Common.get_kernel("FORCE_kernel.cu", block_width, block_height)
self.force_kernel = self.force_module.get_function("FORCEKernel")
self.force_kernel.prepare("iiffffPiPiPiPiPiPi")
#Create data by uploading to device
ghost_cells_x = 1
ghost_cells_y = 1
self.cl_data = Common.SWEDataArkawaA(self.cl_ctx, nx, ny, ghost_cells_x, ghost_cells_y, h0, hu0, hv0)
self.data = Common.SWEDataArakawaA(nx, ny, \
ghost_cells_x, ghost_cells_y, \
h0, hu0, hv0, \
stream=self.stream)
#Save input parameters
#Notice that we need to specify them in the correct dataformat for the
@ -86,7 +92,7 @@ class FORCE:
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]) \
@ -109,20 +115,20 @@ class FORCE:
if (local_dt <= 0.0):
break
self.kernel.swe_2D(self.cl_queue, self.global_size, self.local_size, \
self.force_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.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.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.t += local_dt
self.cl_data.swap()
self.data.swap()
return self.t
@ -131,5 +137,5 @@ class FORCE:
def download(self):
return self.cl_data.download(self.cl_queue)
return self.data.download(self.stream)

70
SWESimulators/FORCE_kernel.opencl → SWESimulators/FORCE_kernel.cu
View File

@ -19,14 +19,15 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "common.opencl"
#include "common.cu"
/**
* Computes the flux along the x axis for all faces
*/
void computeFluxF(__local float Q[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+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
@ -40,12 +41,12 @@ void computeFluxF(__local float Q[3][block_height+2][block_width+2],
const int k = i;
// Q at interface from the right and left
const float3 Qp = (float3)(Q[0][l][k+1],
Q[1][l][k+1],
Q[2][l][k+1]);
const float3 Qm = (float3)(Q[0][l][k],
Q[1][l][k],
Q[2][l][k]);
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_);
@ -54,15 +55,16 @@ void computeFluxF(__local float Q[3][block_height+2][block_width+2],
F[2][j][i] = flux.z;
}
}
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
}
/**
* Computes the flux along the y axis for all faces
*/
void computeFluxG(__local float Q[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+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);
@ -76,12 +78,12 @@ void computeFluxG(__local float Q[3][block_height+2][block_width+2],
// Q at interface from the right and left
// Note that we swap hu and hv
const float3 Qp = (float3)(Q[0][l+1][k],
Q[2][l+1][k],
Q[1][l+1][k]);
const float3 Qm = (float3)(Q[0][l][k],
Q[2][l][k],
Q[1][l][k]);
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
@ -91,24 +93,24 @@ void computeFluxG(__local float Q[3][block_height+2][block_width+2],
G[2][j][i] = flux.y;
}
}
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
}
__kernel void swe_2D(
__global__ void FORCEKernel(
int nx_, int ny_,
float dx_, float dy_, float dt_,
float g_,
//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_) {
//Index of thread within block
const int tx = get_local_id(0);
@ -122,8 +124,8 @@ __kernel void swe_2D(
const int ti = get_global_id(0) + 1; //Skip global ghost cells, i.e., +1
const int tj = get_global_id(1) + 1;
__local float Q[3][block_height+2][block_width+2];
__local float F[3][block_height+1][block_width+1];
__shared__ float Q[3][block_height+2][block_width+2];
__shared__ float F[3][block_height+1][block_width+1];
//Read into shared memory
@ -131,7 +133,7 @@ __kernel void swe_2D(
hu0_ptr_, hu0_pitch_,
hv0_ptr_, hv0_pitch_,
Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Save our input variables
@ -142,23 +144,21 @@ __kernel void swe_2D(
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Compute flux along x, and evolve
computeFluxF(Q, F, g_, dx_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
evolveF1(Q, F, nx_, ny_, dx_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Set boundary conditions
noFlowBoundary1(Q, nx_, ny_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Compute flux along y, and evolve
computeFluxG(Q, F, g_, dy_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
evolveG1(Q, F, nx_, ny_, dy_, dt_);
barrier(CLK_LOCAL_MEM_FENCE);
__syncthreads();
//Write to main memory
writeBlock1(h1_ptr_, h1_pitch_,