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

Fixed order again

Browse Source
This commit is contained in:
André R. Brodtkorb 2018-11-15 16:47:13 +01:00
parent dcb849b705
commit 7592ad5b9f
22 changed files with 758 additions and 619 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
GPUSimulators/Autotuner.py
View File

@ -29,7 +29,7 @@ from socket import gethostname
import pycuda.driver as cuda import pycuda.driver as cuda
from GPUSimulators import Common, Simulator from GPUSimulators import Common, Simulator, CudaContext
class Autotuner: class Autotuner:
def __init__(self, def __init__(self,

28
GPUSimulators/CudaContext.py
View File

@ -48,7 +48,7 @@ class CudaContext(object):
self.blocking = blocking self.blocking = blocking
self.use_cache = use_cache self.use_cache = use_cache
self.logger = logging.getLogger(__name__) self.logger = logging.getLogger(__name__)
self.kernels = {} self.modules = {}
self.module_path = os.path.dirname(os.path.realpath(__file__)) self.module_path = os.path.dirname(os.path.realpath(__file__))
@ -164,12 +164,12 @@ class CudaContext(object):
break break
return kernel_hasher.hexdigest() return kernel_hasher.hexdigest()
""" """
Reads a text file and creates an OpenCL kernel from that Reads a text file and creates an OpenCL kernel from that
""" """
def get_prepared_kernel(self, kernel_filename, kernel_function_name, \ def get_module(self, kernel_filename,
prepared_call_args, \
include_dirs=[], \ include_dirs=[], \
defines={}, \ defines={}, \
compile_args={'no_extern_c', True}, jit_compile_args={}): compile_args={'no_extern_c', True}, jit_compile_args={}):
@ -206,9 +206,9 @@ class CudaContext(object):
cached_kernel_filename = os.path.join(self.cache_path, kernel_hash) cached_kernel_filename = os.path.join(self.cache_path, kernel_hash)
# If we have the kernel in our hashmap, return it # If we have the kernel in our hashmap, return it
if (kernel_hash in self.kernels.keys()): if (kernel_hash in self.modules.keys()):
self.logger.debug("Found kernel %s cached in hashmap (%s)", kernel_filename, kernel_hash) self.logger.debug("Found kernel %s cached in hashmap (%s)", kernel_filename, kernel_hash)
return self.kernels[kernel_hash] return self.modules[kernel_hash]
# If we have it on disk, return it # If we have it on disk, return it
elif (self.use_cache and os.path.isfile(cached_kernel_filename)): elif (self.use_cache and os.path.isfile(cached_kernel_filename)):
@ -218,10 +218,8 @@ class CudaContext(object):
file_str = file.read() file_str = file.read()
module = cuda.module_from_buffer(file_str, message_handler=cuda_compile_message_handler, **jit_compile_args) module = cuda.module_from_buffer(file_str, message_handler=cuda_compile_message_handler, **jit_compile_args)
kernel = module.get_function(kernel_function_name) self.modules[kernel_hash] = module
kernel.prepare(prepared_call_args) return module
self.kernels[kernel_hash] = kernel
return kernel
# Otherwise, compile it from source # Otherwise, compile it from source
else: else:
@ -250,19 +248,15 @@ class CudaContext(object):
with io.open(cached_kernel_filename, "wb") as file: with io.open(cached_kernel_filename, "wb") as file:
file.write(cubin) file.write(cubin)
kernel = module.get_function(kernel_function_name) self.modules[kernel_hash] = module
kernel.prepare(prepared_call_args) return module
self.kernels[kernel_hash] = kernel
return kernel
""" """
Clears the kernel cache (useful for debugging & development) Clears the kernel cache (useful for debugging & development)
""" """
def clear_kernel_cache(self): def clear_kernel_cache(self):
self.logger.debug("Clearing cache") self.logger.debug("Clearing cache")
self.kernels = {} self.modules = {}
gc.collect() gc.collect()
""" """

100
GPUSimulators/EE2D_KP07_dimsplit.py
View File

@ -24,6 +24,8 @@ from GPUSimulators import Simulator, Common
from GPUSimulators.Simulator import BaseSimulator, BoundaryCondition from GPUSimulators.Simulator import BaseSimulator, BoundaryCondition
import numpy as np import numpy as np
from pycuda import gpuarray
@ -52,80 +54,81 @@ class EE2D_KP07_dimsplit (BaseSimulator):
gamma: Gas constant gamma: Gas constant
p: pressure p: pressure
""" """
def __init__(self, \ def __init__(self,
context, \ context,
rho, rho_u, rho_v, E, \ rho, rho_u, rho_v, E,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
gamma, \ gamma,
theta=1.3, \ theta=1.3,
order=2, \ cfl_scale=0.25*0.9,
boundary_conditions=BoundaryCondition(), \ boundary_conditions=BoundaryCondition(),
block_width=16, block_height=8): block_width=16, block_height=8):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context, \
nx, ny, \ nx, ny, \
dx, dy, dt, \ dx, dy, 2*dt, \
block_width, block_height) block_width, block_height)
self.g = np.float32(g) self.g = np.float32(g)
self.gamma = np.float32(gamma) self.gamma = np.float32(gamma)
self.theta = np.float32(theta) self.theta = np.float32(theta)
self.order = np.int32(order) self.cfl_scale = cfl_scale
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/EE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \ module = context.get_module("cuda/EE2D_KP07_dimsplit.cu",
"iiffffffiiPiPiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("KP07DimsplitKernel")
self.kernel.prepare("iiffffffiiPiPiPiPiPiPiPiPiP")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[rho, rho_u, rho_v, E]) [rho, rho_u, rho_v, E])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[None, None, None, None]) [None, None, None, None])
self.cfl_data = gpuarray.GPUArray(self.grid_size, dtype=np.float32)
self.cfl_data.fill(self.dt, stream=self.stream)
def step(self, dt): def step(self, dt):
if (self.order == 1): self.substepDimsplit(0.5*dt, 0)
self.substepDimsplit(dt, substep=(self.nt % 2)) self.substepDimsplit(0.5*dt, 1)
elif (self.order == 2):
self.substepDimsplit(dt, substep=0)
self.substepDimsplit(dt, substep=1)
else:
raise(NotImplementedError("Order {:d} is not implemented".format(self.order)))
self.t += dt self.t += dt
self.nt += 1 self.nt += 2
def substepDimsplit(self, dt, substep): def substepDimsplit(self, dt, substep):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.gamma, \ self.gamma,
self.theta, \ self.theta,
Simulator.stepOrderToCodedInt(step=substep, order=self.order), \ substep,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u0[3].data.gpudata, self.u0[3].data.strides[0], \ self.u0[3].data.gpudata, self.u0[3].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0], \ self.u1[2].data.gpudata, self.u1[2].data.strides[0],
self.u1[3].data.gpudata, self.u1[3].data.strides[0]) self.u1[3].data.gpudata, self.u1[3].data.strides[0],
self.cfl_data.gpudata)
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0
def download(self): def download(self):
@ -134,4 +137,7 @@ class EE2D_KP07_dimsplit (BaseSimulator):
def check(self): def check(self):
self.u0.check() self.u0.check()
self.u1.check() self.u1.check()
pass
def computeDt(self):
max_dt = gpuarray.min(self.cfl_data, stream=self.stream).get();
return max_dt*self.cfl_scale

68
GPUSimulators/FORCE.py
View File

@ -52,61 +52,67 @@ class FORCE (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
boundary_conditions=BoundaryCondition(), \ boundary_conditions=BoundaryCondition(),
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
block_width, block_height); block_width, block_height);
self.g = np.float32(g) self.g = np.float32(g)
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_FORCE.cu", "FORCEKernel", \ module = context.get_module("cuda/SWE2D_FORCE.cu",
"iiffffiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("FORCEKernel")
self.kernel.prepare("iiffffiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
1, 1, \ 1, 1,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
1, 1, \ 1, 1,
[None, None, None]) [None, None, None])
def step(self, dt): def step(self, dt):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0
self.t += dt self.t += dt
self.nt += 1 self.nt += 1
def download(self): def download(self):
return self.u0.download(self.stream) return self.u0.download(self.stream)
def check(self):
self.u0.check()
self.u1.check()

61
GPUSimulators/HLL.py
View File

@ -47,57 +47,58 @@ class HLL (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
boundary_conditions=BoundaryCondition(), \ boundary_conditions=BoundaryCondition(),
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
block_width, block_height); block_width, block_height);
self.g = np.float32(g) self.g = np.float32(g)
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_HLL.cu", "HLLKernel", \ module = context.get_module("cuda/SWE2D_HLL.cu",
"iiffffiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("HLLKernel")
self.kernel.prepare("iiffffiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
1, 1, \ 1, 1,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
1, 1, \ 1, 1,
[None, None, None]) [None, None, None])
def step(self, dt): def step(self, dt):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0
self.t += dt self.t += dt

84
GPUSimulators/HLL2.py
View File

@ -49,79 +49,69 @@ class HLL2 (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
theta=1.8, \ theta=1.8,
order=2, \ boundary_conditions=BoundaryCondition(),
boundary_conditions=BoundaryCondition(), \
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt*2,
block_width, block_height); block_width, block_height);
self.g = np.float32(g) self.g = np.float32(g)
self.theta = np.float32(theta) self.theta = np.float32(theta)
self.order = np.int32(order)
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#This kernel is dimensionally split, and therefore only second order every other
#dimsplit timestep. Therefore, step always runs two substeps
self.dt = 2*self.dt
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_HLL2.cu", "HLL2Kernel", \ module = context.get_module("cuda/SWE2D_HLL2.cu",
"iifffffiiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("HLL2Kernel")
self.kernel.prepare("iifffffiiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[None, None, None]) [None, None, None])
def step(self, dt): def step(self, dt):
if (self.order == 1): self.substepDimsplit(dt*0.5, 0)
self.substepDimsplit(0.5*dt, 0) self.substepDimsplit(dt*0.5, 1)
self.substepDimsplit(0.5*dt, 1)
elif (self.order == 2):
self.substepDimsplit(0.5*dt, 0)
self.substepDimsplit(0.5*dt, 1)
else:
raise(NotImplementedError("Order {:d} is not implemented".format(self.order)))
self.t += dt self.t += dt
self.nt += 2 self.nt += 2
def substepDimsplit(self, dt, substep): def substepDimsplit(self, dt, substep):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.theta, \ self.theta,
Simulator.stepOrderToCodedInt(step=substep, order=self.order), \ substep,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0

72
GPUSimulators/KP07.py
View File

@ -50,21 +50,21 @@ class KP07 (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
theta=1.3, \ theta=1.3,
order=2, \ order=2,
boundary_conditions=BoundaryCondition(), \ boundary_conditions=BoundaryCondition(),
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
block_width, block_height); block_width, block_height);
self.g = np.float32(g) self.g = np.float32(g)
self.theta = np.float32(theta) self.theta = np.float32(theta)
@ -72,26 +72,27 @@ class KP07 (Simulator.BaseSimulator):
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_KP07.cu", "KP07Kernel", \ module = context.get_module("cuda/SWE2D_KP07.cu",
"iifffffiiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("KP07Kernel")
self.kernel.prepare("iifffffiiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[None, None, None]) [None, None, None])
@ -108,20 +109,21 @@ class KP07 (Simulator.BaseSimulator):
def substepRK(self, dt, substep): def substepRK(self, dt, substep):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.theta, \ self.theta,
Simulator.stepOrderToCodedInt(step=substep, order=self.order), \ Simulator.stepOrderToCodedInt(step=substep, order=self.order),
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0
def download(self): def download(self):
return self.u0.download(self.stream) return self.u0.download(self.stream)

90
GPUSimulators/KP07_dimsplit.py
View File

@ -50,76 +50,78 @@ class KP07_dimsplit (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
theta=1.3, \ theta=1.3,
order=2, \ boundary_conditions=BoundaryCondition(),
boundary_conditions=BoundaryCondition(), \
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt*2,
block_width, block_height) block_width, block_height)
self.gc_x = 2
self.gc_y = 2
self.g = np.float32(g) self.g = np.float32(g)
self.theta = np.float32(theta) self.theta = np.float32(theta)
self.order = np.int32(order)
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \ module = context.get_module("cuda/SWE2D_KP07_dimsplit.cu",
"iifffffiiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("KP07DimsplitKernel")
self.kernel.prepare("iifffffiiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ self.gc_x, self.gc_y,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ self.gc_x, self.gc_y,
[None, None, None]) [None, None, None])
def step(self, dt): def step(self, dt):
if (self.order == 1): self.substepDimsplit(dt*0.5, 0)
self.substepDimsplit(dt, substep=(self.nt % 2)) self.substepDimsplit(dt*0.5, 1)
elif (self.order == 2):
self.substepDimsplit(dt, substep=0)
self.substepDimsplit(dt, substep=1)
else:
raise(NotImplementedError("Order {:d} is not implemented".format(self.order)))
self.t += dt self.t += dt
self.nt += 1 self.nt += 2
def substepDimsplit(self, dt, substep): def substepDimsplit(self, dt, substep):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.theta, \ self.theta,
Simulator.stepOrderToCodedInt(step=substep, order=self.order), \ substep,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0
def download(self): def download(self):
return self.u0.download(self.stream) return self.u0.download(self.stream)
def check(self):
self.u0.check()
self.u1.check()

59
GPUSimulators/LxF.py
View File

@ -48,57 +48,58 @@ class LxF (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
boundary_conditions=BoundaryCondition(), boundary_conditions=BoundaryCondition(),
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
block_width, block_height); block_width, block_height);
self.g = np.float32(g) self.g = np.float32(g)
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
# Get kernels # Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_LxF.cu", "LxFKernel", \ module = context.get_module("cuda/SWE2D_LxF.cu",
"iiffffiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("LxFKernel")
self.kernel.prepare("iiffffiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
1, 1, \ 1, 1,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
1, 1, \ 1, 1,
[None, None, None]) [None, None, None])
def step(self, dt): def step(self, dt):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0
self.t += dt self.t += dt

60
GPUSimulators/Simulator.py
View File

@ -50,9 +50,9 @@ class BoundaryCondition(object):
Neumann = 1, Neumann = 1,
Periodic = 2, Periodic = 2,
Reflective = 3 Reflective = 3
def __init__(self, types={ \ def __init__(self, types={ \
'north': Type.Reflective, \ 'north': Type.Reflective, \
'south': Type.Reflective, \ 'south': Type.Reflective, \
@ -85,13 +85,13 @@ class BoundaryCondition(object):
bc = bc | (self.north & 0x0000000F) << 24 bc = bc | (self.north & 0x0000000F) << 24
bc = bc | (self.south & 0x0000000F) << 16 bc = bc | (self.south & 0x0000000F) << 16
bc = bc | (self.east & 0x0000000F) << 8 bc = bc | (self.east & 0x0000000F) << 8
bc = bc | (self.west & 0x0000000F) bc = bc | (self.west & 0x0000000F)
#for t in types: #for t in types:
# print("{0:s}, {1:d}, {1:032b}, {1:08b}".format(t, types[t])) # print("{0:s}, {1:d}, {1:032b}, {1:08b}".format(t, types[t]))
#print("bc: {0:032b}".format(bc)) #print("bc: {0:032b}".format(bc))
return np.int32(bc) return np.int32(bc)
@ -101,10 +101,10 @@ class BoundaryCondition(object):
class BaseSimulator(object): class BaseSimulator(object):
def __init__(self, \ def __init__(self,
context, \ context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
block_width, block_height): block_width, block_height):
""" """
Initialization routine Initialization routine
@ -141,9 +141,9 @@ class BaseSimulator(object):
#Compute kernel launch parameters #Compute kernel launch parameters
self.block_size = (block_width, block_height, 1) self.block_size = (block_width, block_height, 1)
self.grid_size = ( \ self.grid_size = (
int(np.ceil(self.nx / float(self.block_size[0]))), \ int(np.ceil(self.nx / float(self.block_size[0]))),
int(np.ceil(self.ny / float(self.block_size[1]))) \ int(np.ceil(self.ny / float(self.block_size[1])))
) )
#Create a CUDA stream #Create a CUDA stream
@ -158,43 +158,42 @@ class BaseSimulator(object):
return "{:s} [{:d}x{:d}]".format(self.__class__.__name__, self.nx, self.ny) return "{:s} [{:d}x{:d}]".format(self.__class__.__name__, self.nx, self.ny)
def simulate(self, t_end): def simulate(self, t):
""" """
Function which simulates t_end seconds using the step function Function which simulates t_end seconds using the step function
Requires that the step() function is implemented in the subclasses Requires that the step() function is implemented in the subclasses
""" """
# Compute number of timesteps to perform
n = int(t_end / self.dt + 1)
printer = Common.ProgressPrinter(n) printer = Common.ProgressPrinter(t)
t_end = self.simTime() + t
while(self.simTime() < t_end):
if (self.simSteps() % 100 == 0):
self.dt = self.computeDt()
for i in range(0, n):
# Compute timestep for "this" iteration (i.e., shorten last timestep) # Compute timestep for "this" iteration (i.e., shorten last timestep)
local_dt = np.float32(min(self.dt, t_end-i*self.dt)) local_dt = np.float32(min(self.dt, t_end-self.simTime()))
# Stop if end reached (should not happen) # Stop if end reached (should not happen)
if (local_dt <= 0.0): if (local_dt <= 0.0):
self.logger.warning("Timestep size {:d} is less than or equal to zero!".format(self.nt + i)) self.logger.warning("Timestep size {:d} is less than or equal to zero!".format(self.simSteps()))
break break
# Step forward in time # Step forward in time
self.step(local_dt) self.step(local_dt)
#Print info #Print info
print_string = printer.getPrintString(i) print_string = printer.getPrintString(t_end - self.simTime())
if (print_string): if (print_string):
self.logger.info("%s (Euler): %s", self, print_string) self.logger.info("%s: %s", self, print_string)
try: try:
self.check() self.check()
except AssertionError as e: except AssertionError as e:
e.args += ("Step={:d}, time={:f}".format(self.simSteps(), self.simTime())) e.args += ("Step={:d}, time={:f}".format(self.simSteps(), self.simTime()))
raise raise
#self.logger.info("%s simulated %f seconds to %f with %d steps (Euler)", self, t_end, self.t, n)
return self.t, n
def step(self, dt): def step(self, dt):
""" """
Function which performs one single timestep of size dt Function which performs one single timestep of size dt
@ -208,7 +207,8 @@ class BaseSimulator(object):
self.stream.synchronize() self.stream.synchronize()
def check(self): def check(self):
raise(NotImplementedError("Needs to be implemented in subclass")) self.logger.warning("check() is not implemented - please implement")
#raise(NotImplementedError("Needs to be implemented in subclass"))
def simTime(self): def simTime(self):
return self.t return self.t
@ -216,6 +216,8 @@ class BaseSimulator(object):
def simSteps(self): def simSteps(self):
return self.nt return self.nt
def computeDt(self):
raise(NotImplementedError("Needs to be implemented in subclass"))

76
GPUSimulators/WAF.py
View File

@ -46,71 +46,65 @@ class WAF (Simulator.BaseSimulator):
dt: Size of each timestep (90 s) dt: Size of each timestep (90 s)
g: Gravitational accelleration (9.81 m/s^2) g: Gravitational accelleration (9.81 m/s^2)
""" """
def __init__(self, \ def __init__(self,
context, \ context,
h0, hu0, hv0, \ h0, hu0, hv0,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt,
g, \ g,
order=2, \ boundary_conditions=BoundaryCondition(),
boundary_conditions=BoundaryCondition(), \
block_width=16, block_height=16): block_width=16, block_height=16):
# Call super constructor # Call super constructor
super().__init__(context, \ super().__init__(context,
nx, ny, \ nx, ny,
dx, dy, dt, \ dx, dy, dt*2,
block_width, block_height); block_width, block_height);
self.g = np.float32(g) self.g = np.float32(g)
self.order = np.int32(order)
self.boundary_conditions = boundary_conditions.asCodedInt() self.boundary_conditions = boundary_conditions.asCodedInt()
#Get kernels #Get kernels
self.kernel = context.get_prepared_kernel("cuda/SWE2D_WAF.cu", "WAFKernel", \ module = context.get_module("cuda/SWE2D_WAF.cu",
"iiffffiiPiPiPiPiPiPi", \
defines={ defines={
'BLOCK_WIDTH': self.block_size[0], 'BLOCK_WIDTH': self.block_size[0],
'BLOCK_HEIGHT': self.block_size[1] 'BLOCK_HEIGHT': self.block_size[1]
}, \ },
compile_args={ compile_args={
'no_extern_c': True, 'no_extern_c': True,
'options': ["--use_fast_math"], 'options': ["--use_fast_math"],
}, \ },
jit_compile_args={}) jit_compile_args={})
self.kernel = module.get_function("WAFKernel")
self.kernel.prepare("iiffffiiPiPiPiPiPiPi")
#Create data by uploading to device #Create data by uploading to device
self.u0 = Common.ArakawaA2D(self.stream, \ self.u0 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[h0, hu0, hv0]) [h0, hu0, hv0])
self.u1 = Common.ArakawaA2D(self.stream, \ self.u1 = Common.ArakawaA2D(self.stream,
nx, ny, \ nx, ny,
2, 2, \ 2, 2,
[None, None, None]) [None, None, None])
def step(self, dt): def step(self, dt):
if (self.order == 1): self.substepDimsplit(dt*0.5, substep=0)
self.substepDimsplit(dt, substep=(self.nt % 2)) self.substepDimsplit(dt*0.5, substep=1)
elif (self.order == 2):
self.substepDimsplit(dt, substep=0)
self.substepDimsplit(dt, substep=1)
else:
raise(NotImplementedError("Order {:d} is not implemented".format(self.order)))
self.t += dt self.t += dt
self.nt += 1 self.nt += 2
def substepDimsplit(self, dt, substep): def substepDimsplit(self, dt, substep):
self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream,
self.nx, self.ny, \ self.nx, self.ny,
self.dx, self.dy, dt, \ self.dx, self.dy, dt,
self.g, \ self.g,
Simulator.stepOrderToCodedInt(step=substep, order=self.order), \ substep,
self.boundary_conditions, \ self.boundary_conditions,
self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[0].data.gpudata, self.u0[0].data.strides[0],
self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0],
self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0],
self.u1[0].data.gpudata, self.u1[0].data.strides[0], \ self.u1[0].data.gpudata, self.u1[0].data.strides[0],
self.u1[1].data.gpudata, self.u1[1].data.strides[0], \ self.u1[1].data.gpudata, self.u1[1].data.strides[0],
self.u1[2].data.gpudata, self.u1[2].data.strides[0]) self.u1[2].data.gpudata, self.u1[2].data.strides[0])
self.u0, self.u1 = self.u1, self.u0 self.u0, self.u1 = self.u1, self.u0

80
GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
View File

@ -58,8 +58,8 @@ void computeFluxF(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_)); const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
// Compute flux based on prediction // Compute flux based on prediction
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_); //const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
//const float4 flux = HLL_flux(Q_l_bar, Q_r_bar, gamma_); const float4 flux = HLL_flux(Q_l_bar, Q_r_bar, gamma_);
//Write to shared memory //Write to shared memory
F[0][j][i] = flux.x; F[0][j][i] = flux.x;
@ -131,7 +131,7 @@ __global__ void KP07DimsplitKernel(
float theta_, float theta_,
int step_order_, int step_,
int boundary_conditions_, int boundary_conditions_,
//Input h^n //Input h^n
@ -144,51 +144,49 @@ __global__ void KP07DimsplitKernel(
float* rho1_ptr_, int rho1_pitch_, float* rho1_ptr_, int rho1_pitch_,
float* rho_u1_ptr_, int rho_u1_pitch_, float* rho_u1_ptr_, int rho_u1_pitch_,
float* rho_v1_ptr_, int rho_v1_pitch_, float* rho_v1_ptr_, int rho_v1_pitch_,
float* E1_ptr_, int E1_pitch_) { float* E1_ptr_, int E1_pitch_,
//Output CFL
float* cfl_) {
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 2; const unsigned int gc_x = 2;
const unsigned int gc_y = 2;
const unsigned int vars = 4; const unsigned int vars = 4;
//Shared memory variables //Shared memory variables
__shared__ float Q[4][h+4][w+4]; __shared__ float Q[4][h+2*gc_y][w+2*gc_x];
__shared__ float Qx[4][h+4][w+4]; __shared__ float Qx[4][h+2*gc_y][w+2*gc_x];
__shared__ float F[4][h+4][w+4]; __shared__ float F[4][h+2*gc_y][w+2*gc_x];
//Read into shared memory //Read into shared memory
readBlock<w, h, gc, 1, 1>( rho0_ptr_, rho0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( rho0_ptr_, rho0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(rho_u0_ptr_, rho_u0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(rho_u0_ptr_, rho_u0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(rho_v0_ptr_, rho_v0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(rho_v0_ptr_, rho_v0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, 1>( E0_ptr_, E0_pitch_, Q[3], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( E0_ptr_, E0_pitch_, Q[3], nx_, ny_, boundary_conditions_);
__syncthreads();
//Step 0 => evolve x first, then y //Step 0 => evolve x first, then y
if (getStep(step_order_) == 0) { if (step_ == 0) {
//Compute fluxes along the x axis and evolve //Compute fluxes along the x axis and evolve
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeX<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxF(Q, Qx, F, gamma_, dx_, dt_); computeFluxF(Q, Qx, F, gamma_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
evolveF<w, h, gc, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the y axis and evolve //Compute fluxes along the y axis and evolve
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeY<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxG(Q, Qx, F, gamma_, dy_, dt_); computeFluxG(Q, Qx, F, gamma_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
//Gravity source term //Gravity source term
if (g_ > 0.0f) { if (g_ > 0.0f) {
const int i = threadIdx.x + gc; const int i = threadIdx.x + gc_x;
const int j = threadIdx.y + gc; const int j = threadIdx.y + gc_y;
const float rho_v = Q[2][j][i]; const float rho_v = Q[2][j][i];
Q[2][j][i] -= g_*Q[0][j][i]*dt_; Q[2][j][i] -= g_*Q[0][j][i]*dt_;
Q[3][j][i] -= g_*rho_v*dt_; Q[3][j][i] -= g_*rho_v*dt_;
@ -198,29 +196,25 @@ __global__ void KP07DimsplitKernel(
//Step 1 => evolve y first, then x //Step 1 => evolve y first, then x
else { else {
//Compute fluxes along the y axis and evolve //Compute fluxes along the y axis and evolve
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeY<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxG(Q, Qx, F, gamma_, dy_, dt_); computeFluxG(Q, Qx, F, gamma_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the x axis and evolve //Compute fluxes along the x axis and evolve
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeX<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxF(Q, Qx, F, gamma_, dx_, dt_); computeFluxF(Q, Qx, F, gamma_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
evolveF<w, h, gc, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Gravity source term //Gravity source term
if (g_ > 0.0f) { if (g_ > 0.0f) {
const int i = threadIdx.x + gc; const int i = threadIdx.x + gc_x;
const int j = threadIdx.y + gc; const int j = threadIdx.y + gc_y;
const float rho_v = Q[2][j][i]; const float rho_v = Q[2][j][i];
Q[2][j][i] -= g_*Q[0][j][i]*dt_; Q[2][j][i] -= g_*Q[0][j][i]*dt_;
Q[3][j][i] -= g_*rho_v*dt_; Q[3][j][i] -= g_*rho_v*dt_;
@ -230,12 +224,16 @@ __global__ void KP07DimsplitKernel(
// Write to main memory for all internal cells // Write to main memory for all internal cells
const int step = getStep(step_order_); writeBlock<w, h, gc_x, gc_y>( rho1_ptr_, rho1_pitch_, Q[0], nx_, ny_, 0, 1);
const int order = getOrder(step_order_); writeBlock<w, h, gc_x, gc_y>(rho_u1_ptr_, rho_u1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, gc>( rho1_ptr_, rho1_pitch_, Q[0], nx_, ny_, step, order); writeBlock<w, h, gc_x, gc_y>(rho_v1_ptr_, rho_v1_pitch_, Q[2], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(rho_u1_ptr_, rho_u1_pitch_, Q[1], nx_, ny_, step, order); writeBlock<w, h, gc_x, gc_y>( E1_ptr_, E1_pitch_, Q[3], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(rho_v1_ptr_, rho_v1_pitch_, Q[2], nx_, ny_, step, order);
writeBlock<w, h, gc>( E1_ptr_, E1_pitch_, Q[3], nx_, ny_, step, order); //Compute the CFL for this block
if (cfl_ != NULL) {
writeCfl<w, h, gc_x, gc_y, vars>(Q, F[0], nx_, ny_, dx_, dy_, gamma_, cfl_);
}
} }
} // extern "C" } // extern "C"

54
GPUSimulators/cuda/EulerCommon.h
View File

@ -23,6 +23,60 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
template<int w, int h, int gc_x, int gc_y, int vars>
__device__ void writeCfl(float Q[vars][h+2*gc_y][w+2*gc_x],
float shmem[h+2*gc_y][w+2*gc_x],
const int nx_, const int ny_,
const float dx_, const float dy_, const float gamma_,
float* output_) {
//Index of thread within block
const int tx = threadIdx.x + gc_x;
const int ty = threadIdx.y + gc_y;
//Index of cell within domain
const int ti = blockDim.x*blockIdx.x + tx;
const int tj = blockDim.y*blockIdx.y + ty;
//Only internal cells
if (ti < nx_+gc_x && tj < ny_+gc_y) {
const float rho = Q[0][ty][tx];
const float u = Q[1][ty][tx] / rho;
const float v = Q[2][ty][tx] / rho;
const float max_u = dx_ / (fabsf(u) + sqrtf(gamma_*rho));
const float max_v = dy_ / (fabsf(v) + sqrtf(gamma_*rho));
shmem[ty][tx] = fminf(max_u, max_v);
}
__syncthreads();
//One row of threads loop over all rows
if (ti < nx_+gc_x && tj < ny_+gc_y) {
if (ty == gc_y) {
float min_val = shmem[ty][tx];
const int max_y = min(h, ny_+gc_y - tj);
for (int j=gc_y; j<max_y+gc_y; j++) {
min_val = fminf(min_val, shmem[j][tx]);
}
shmem[ty][tx] = min_val;
}
}
__syncthreads();
//One thread loops over first row to find global max
if (tx == gc_x && ty == gc_y) {
float min_val = shmem[ty][tx];
const int max_x = min(w, nx_+gc_x - ti);
for (int i=gc_x; i<max_x+gc_x; ++i) {
min_val = fminf(min_val, shmem[ty][i]);
}
const int idx = gridDim.x*blockIdx.y + blockIdx.x;
output_[idx] = min_val;
}
}

23
GPUSimulators/cuda/SWE2D_FORCE.cu
View File

@ -101,34 +101,35 @@ __global__ void FORCEKernel(
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 1; const unsigned int gc_x = 1;
const unsigned int gc_y = 1;
const unsigned int vars = 3; const unsigned int vars = 3;
__shared__ float Q[3][h+2][w+2]; __shared__ float Q[vars][h+2*gc_y][w+2*gc_x];
__shared__ float F[3][h+2][w+2]; __shared__ float F[vars][h+2*gc_y][w+2*gc_x];
//Read into shared memory //Read into shared memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
__syncthreads(); __syncthreads();
//Compute flux along x, and evolve //Compute flux along x, and evolve
computeFluxF(Q, F, g_, dx_, dt_); computeFluxF(Q, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Compute flux along y, and evolve //Compute flux along y, and evolve
computeFluxG(Q, F, g_, dy_, dt_); computeFluxG(Q, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
//Write to main memory //Write to main memory
writeBlock<w, h, gc>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1);
} }
} // extern "C" } // extern "C"

23
GPUSimulators/cuda/SWE2D_HLL.cu
View File

@ -117,36 +117,37 @@ __global__ void HLLKernel(
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 1; const unsigned int gc_x = 1;
const unsigned int gc_y = 1;
const unsigned int vars = 3; const unsigned int vars = 3;
//Shared memory variables //Shared memory variables
__shared__ float Q[3][h+2][w+2]; __shared__ float Q[vars][h+2*gc_y][w+2*gc_x];
__shared__ float F[3][h+2][w+2]; __shared__ float F[vars][h+2*gc_y][w+2*gc_x];
//Read into shared memory //Read into shared memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
//Compute F flux //Compute F flux
computeFluxF(Q, F, g_); computeFluxF(Q, F, g_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Compute G flux //Compute G flux
computeFluxG(Q, F, g_); computeFluxG(Q, F, g_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
// Write to main memory for all internal cells // Write to main memory for all internal cells
writeBlock<w, h, gc>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1);
} }
} // extern "C" } // extern "C"

37
GPUSimulators/cuda/SWE2D_HLL2.cu
View File

@ -130,7 +130,7 @@ __global__ void HLL2Kernel(
float theta_, float theta_,
int step_order_, int step_,
int boundary_conditions_, int boundary_conditions_,
//Input h^n //Input h^n
@ -145,7 +145,8 @@ __global__ void HLL2Kernel(
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 2; const unsigned int gc_x = 2;
const unsigned int gc_y = 2;
const unsigned int vars = 3; const unsigned int vars = 3;
//Shared memory variables //Shared memory variables
@ -154,44 +155,44 @@ __global__ void HLL2Kernel(
__shared__ float F[3][h+4][w+4]; __shared__ float F[3][h+4][w+4];
//Read into shared memory //Read into shared memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
//Step 0 => evolve x first, then y //Step 0 => evolve x first, then y
if (getStep(step_order_) == 0) { if (step_ == 0) {
//Compute fluxes along the x axis and evolve //Compute fluxes along the x axis and evolve
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeX<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_); computeFluxF(Q, Qx, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the y axis and evolve //Compute fluxes along the y axis and evolve
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeY<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_); computeFluxG(Q, Qx, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
} }
//Step 1 => evolve y first, then x //Step 1 => evolve y first, then x
else { else {
//Compute fluxes along the y axis and evolve //Compute fluxes along the y axis and evolve
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeY<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_); computeFluxG(Q, Qx, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the x axis and evolve //Compute fluxes along the x axis and evolve
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeX<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_); computeFluxF(Q, Qx, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
} }
@ -199,11 +200,9 @@ __global__ void HLL2Kernel(
// Write to main memory for all internal cells // Write to main memory for all internal cells
const int step = getStep(step_order_); writeBlock<w, h, gc_x, gc_y>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1);
const int order = getOrder(step_order_); writeBlock<w, h, gc_x, gc_y>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, 2>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, step, order); writeBlock<w, h, gc_x, gc_y>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1);
writeBlock<w, h, 2>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, step, order);
writeBlock<w, h, 2>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, step, order);
} }
} // extern "C" } // extern "C"

9
GPUSimulators/cuda/SWE2D_KP07.cu
View File

@ -155,7 +155,8 @@ __global__ void KP07Kernel(
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 2; const unsigned int gc_x = 2;
const unsigned int gc_y = 2;
const unsigned int vars = 3; const unsigned int vars = 3;
//Index of thread within block //Index of thread within block
@ -176,9 +177,9 @@ __global__ void KP07Kernel(
//Read into shared memory //Read into shared memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
//Reconstruct slopes along x and axis //Reconstruct slopes along x and axis

103
GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
View File

@ -29,13 +29,14 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "limiters.h" #include "limiters.h"
template <int w, int h, int gc_x, int gc_y>
__device__ __device__
void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], void computeFluxF(float Q[3][h+2*gc_y][w+2*gc_x],
float Qx[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], float Qx[3][h+2*gc_y][w+2*gc_x],
float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], float F[3][h+2*gc_y][w+2*gc_x],
const float g_, const float dx_, const float dt_) { const float g_, const float dx_, const float dt_) {
for (int j=threadIdx.y; j<BLOCK_HEIGHT+4; j+=BLOCK_HEIGHT) { for (int j=threadIdx.y; j<h+2*gc_y; j+=h) {
for (int i=threadIdx.x+1; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) { for (int i=threadIdx.x+1; i<w+2*gc_x-2; i+=w) {
// Reconstruct point values of Q at the left and right hand side // 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 // of the cell for both the left (i) and right (i+1) cell
const float3 Q_rl = make_float3(Q[0][j][i+1] - 0.5f*Qx[0][j][i+1], const float3 Q_rl = make_float3(Q[0][j][i+1] - 0.5f*Qx[0][j][i+1],
@ -67,13 +68,14 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
} }
} }
template <int w, int h, int gc_x, int gc_y>
__device__ __device__
void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], void computeFluxG(float Q[3][h+2*gc_y][w+2*gc_x],
float Qy[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], float Qy[3][h+2*gc_y][w+2*gc_x],
float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], float G[3][h+2*gc_y][w+2*gc_x],
const float g_, const float dy_, const float dt_) { const float g_, const float dy_, const float dt_) {
for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) { for (int j=threadIdx.y+1; j<h+2*gc_y-2; j+=h) {
for (int i=threadIdx.x; i<BLOCK_WIDTH+4; i+=BLOCK_WIDTH) { for (int i=threadIdx.x; i<w+2*gc_x; i+=w) {
// Reconstruct point values of Q at the left and right hand side // 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 // of the cell for both the left (i) and right (i+1) cell
//NOte that hu and hv are swapped ("transposing" the domain)! //NOte that hu and hv are swapped ("transposing" the domain)!
@ -114,6 +116,10 @@ void computeFluxG(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
*/ */
extern "C" { extern "C" {
__global__ void KP07DimsplitKernel( __global__ void KP07DimsplitKernel(
int nx_, int ny_, int nx_, int ny_,
float dx_, float dy_, float dt_, float dx_, float dy_, float dt_,
@ -121,7 +127,7 @@ __global__ void KP07DimsplitKernel(
float theta_, float theta_,
int step_order_, int step_,
int boundary_conditions_, int boundary_conditions_,
//Input h^n //Input h^n
@ -133,71 +139,70 @@ __global__ void KP07DimsplitKernel(
float* h1_ptr_, int h1_pitch_, float* h1_ptr_, int h1_pitch_,
float* hu1_ptr_, int hu1_pitch_, float* hu1_ptr_, int hu1_pitch_,
float* hv1_ptr_, int hv1_pitch_) { float* hv1_ptr_, int hv1_pitch_) {
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 2; const unsigned int gc_x = 2;
const unsigned int gc_y = 2;
const unsigned int vars = 3; const unsigned int vars = 3;
//Shared memory variables //Shared memory variables
__shared__ float Q[3][h+4][w+4]; __shared__ float Q[vars][h+2*gc_y][w+2*gc_x];
__shared__ float Qx[3][h+4][w+4]; __shared__ float Qx[vars][h+2*gc_y][w+2*gc_x];
__shared__ float F[3][h+4][w+4]; __shared__ float F[vars][h+2*gc_y][w+2*gc_x];
//Read into shared memory //Read into shared memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
if (step_ == 0) {
//Along X
//Step 0 => evolve x first, then y minmodSlopeX<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
if (getStep(step_order_) == 0) {
//Compute fluxes along the x axis and evolve
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_); computeFluxF<w, h, gc_x, gc_y>(Q, Qx, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the y axis and evolve //Along Y
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeY<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_); computeFluxG<w, h, gc_x, gc_y>(Q, Qx, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
} }
//Step 1 => evolve y first, then x
else { else {
//Compute fluxes along the y axis and evolve //Along Y
minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeY<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxG(Q, Qx, F, g_, dy_, dt_); computeFluxG<w, h, gc_x, gc_y>(Q, Qx, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the x axis and evolve //Along X
minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_); minmodSlopeX<w, h, gc_x, gc_y, vars>(Q, Qx, theta_);
__syncthreads(); __syncthreads();
computeFluxF(Q, Qx, F, g_, dx_, dt_); computeFluxF<w, h, gc_x, gc_y>(Q, Qx, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
} }
// Write to main memory for all internal cells // Write to main memory for all internal cells
const int step = getStep(step_order_); writeBlock<w, h, gc_x, gc_y>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1);
const int order = getOrder(step_order_); writeBlock<w, h, gc_x, gc_y>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, gc>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, step, order); writeBlock<w, h, gc_x, gc_y>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, step, order);
writeBlock<w, h, gc>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, step, order);
} }
} // extern "C" } // extern "C"

22
GPUSimulators/cuda/SWE2D_LxF.cu
View File

@ -118,16 +118,18 @@ void LxFKernel(
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 1; const unsigned int gc_x = 1;
const unsigned int gc_y = 1;
const unsigned int vars = 3;
__shared__ float Q[3][h+2][w+2]; __shared__ float Q[vars][h+2][w+2];
__shared__ float F[3][h ][w+1]; __shared__ float F[vars][h ][w+1];
__shared__ float G[3][h+1][w ]; __shared__ float G[vars][h+1][w ];
//Read from global memory //Read from global memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
//Compute fluxes along the x and y axis //Compute fluxes along the x and y axis
computeFluxF<w, h>(Q, F, g_, dx_, dt_); computeFluxF<w, h>(Q, F, g_, dx_, dt_);
@ -149,9 +151,9 @@ void LxFKernel(
__syncthreads(); __syncthreads();
//Write to main memory //Write to main memory
writeBlock<w, h, gc>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1); writeBlock<w, h, gc_x, gc_y>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1);
} }
} // extern "C" } // extern "C"

29
GPUSimulators/cuda/SWE2D_WAF.cu
View File

@ -105,7 +105,7 @@ __global__ void WAFKernel(
float dx_, float dy_, float dt_, float dx_, float dy_, float dt_,
float g_, float g_,
int step_order_, int step_,
int boundary_conditions_, int boundary_conditions_,
//Input h^n //Input h^n
@ -120,7 +120,8 @@ __global__ void WAFKernel(
const unsigned int w = BLOCK_WIDTH; const unsigned int w = BLOCK_WIDTH;
const unsigned int h = BLOCK_HEIGHT; const unsigned int h = BLOCK_HEIGHT;
const unsigned int gc = 2; const unsigned int gc_x = 2;
const unsigned int gc_y = 2;
const unsigned int vars = 3; const unsigned int vars = 3;
//Shared memory variables //Shared memory variables
@ -130,25 +131,25 @@ __global__ void WAFKernel(
//Read into shared memory Q from global memory //Read into shared memory Q from global memory
readBlock<w, h, gc, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, 1>( h0_ptr_, h0_pitch_, Q[0], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, -1, 1>(hu0_ptr_, hu0_pitch_, Q[1], nx_, ny_, boundary_conditions_);
readBlock<w, h, gc, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_); readBlock<w, h, gc_x, gc_y, 1, -1>(hv0_ptr_, hv0_pitch_, Q[2], nx_, ny_, boundary_conditions_);
__syncthreads(); __syncthreads();
//Step 0 => evolve x first, then y //Step 0 => evolve x first, then y
if (getStep(step_order_) == 0) { if (step_ == 0) {
//Compute fluxes along the x axis and evolve //Compute fluxes along the x axis and evolve
computeFluxF(Q, F, g_, dx_, dt_); computeFluxF(Q, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the y axis and evolve //Compute fluxes along the y axis and evolve
computeFluxG(Q, F, g_, dy_, dt_); computeFluxG(Q, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
} }
//Step 1 => evolve y first, then x //Step 1 => evolve y first, then x
@ -156,24 +157,22 @@ __global__ void WAFKernel(
//Compute fluxes along the y axis and evolve //Compute fluxes along the y axis and evolve
computeFluxG(Q, F, g_, dy_, dt_); computeFluxG(Q, F, g_, dy_, dt_);
__syncthreads(); __syncthreads();
evolveG<w, h, gc, vars>(Q, F, dy_, dt_); evolveG<w, h, gc_x, gc_y, vars>(Q, F, dy_, dt_);
__syncthreads(); __syncthreads();
//Compute fluxes along the x axis and evolve //Compute fluxes along the x axis and evolve
computeFluxF(Q, F, g_, dx_, dt_); computeFluxF(Q, F, g_, dx_, dt_);
__syncthreads(); __syncthreads();
evolveF<w, h, gc, vars>(Q, F, dx_, dt_); evolveF<w, h, gc_x, gc_y, vars>(Q, F, dx_, dt_);
__syncthreads(); __syncthreads();
} }
// Write to main memory for all internal cells // Write to main memory for all internal cells
const int step = getStep(step_order_); writeBlock<w, h, gc_x, gc_y>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, 0, 1);
const int order = getOrder(step_order_); writeBlock<w, h, gc_x, gc_y>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, 0, 1);
writeBlock<w, h, gc>( h1_ptr_, h1_pitch_, Q[0], nx_, ny_, step, order); writeBlock<w, h, gc_x, gc_y>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, 0, 1);
writeBlock<w, h, gc>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_, step, order);
writeBlock<w, h, gc>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_, step, order);
} }
} // extern "C" } // extern "C"

277
GPUSimulators/cuda/common.h
View File

@ -137,19 +137,21 @@ inline __device__ BoundaryCondition getBCWest(int bc_) {
/** /**
* Alter the index l so that it gives periodic boundary conditions when reading * Alter the index l so that it gives periodic boundary conditions when reading
*/ */
template<int ghost_cells> template<int gc_x>
inline __device__ int handlePeriodicBoundaryX(int k, int nx_, int boundary_conditions_) { inline __device__ int handlePeriodicBoundaryX(int k, int nx_, int boundary_conditions_) {
const int gc_pad = 2*ghost_cells; const int gc_pad = gc_x;
//West boundary: add an offset to read from east of domain //West boundary: add an offset to read from east of domain
if ((k < gc_pad) if (gc_x > 0) {
&& getBCWest(boundary_conditions_) == Periodic) { if ((k < gc_pad)
k += (nx_+2*ghost_cells - 2*gc_pad); && getBCWest(boundary_conditions_) == Periodic) {
} k += (nx_+2*gc_x - 2*gc_pad);
//East boundary: subtract an offset to read from west of domain }
else if ((k >= nx_+2*ghost_cells-gc_pad) //East boundary: subtract an offset to read from west of domain
&& getBCEast(boundary_conditions_) == Periodic) { else if ((k >= nx_+2*gc_x-gc_pad)
k -= (nx_+2*ghost_cells - 2*gc_pad); && getBCEast(boundary_conditions_) == Periodic) {
k -= (nx_+2*gc_x - 2*gc_pad);
}
} }
return k; return k;
@ -158,45 +160,49 @@ inline __device__ int handlePeriodicBoundaryX(int k, int nx_, int boundary_condi
/** /**
* Alter the index l so that it gives periodic boundary conditions when reading * Alter the index l so that it gives periodic boundary conditions when reading
*/ */
template<int ghost_cells> template<int gc_y>
inline __device__ int handlePeriodicBoundaryY(int l, int ny_, int boundary_conditions_) { inline __device__ int handlePeriodicBoundaryY(int l, int ny_, int boundary_conditions_) {
const int gc_pad = 2*ghost_cells; const int gc_pad = gc_y;
//South boundary: add an offset to read from north of domain //South boundary: add an offset to read from north of domain
if ((l < gc_pad) if (gc_y > 0) {
&& getBCSouth(boundary_conditions_) == Periodic) { if ((l < gc_pad)
l += (ny_+2*ghost_cells - 2*gc_pad); && getBCSouth(boundary_conditions_) == Periodic) {
} l += (ny_+2*gc_y - 2*gc_pad);
//North boundary: subtract an offset to read from south of domain }
else if ((l >= ny_+2*ghost_cells-gc_pad) //North boundary: subtract an offset to read from south of domain
&& getBCNorth(boundary_conditions_) == Periodic) { else if ((l >= ny_+2*gc_y-gc_pad)
l -= (ny_+2*ghost_cells - 2*gc_pad); && getBCNorth(boundary_conditions_) == Periodic) {
l -= (ny_+2*gc_y - 2*gc_pad);
}
} }
return l; return l;
} }
template<int block_width, int block_height, int ghost_cells, int sign_x, int sign_y> template<int w, int h, int gc_x, int gc_y, int sign_x, int sign_y>
inline __device__ int handleReflectiveBoundary( inline __device__
float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], void handleReflectiveBoundary(
float Q[h+2*gc_y][w+2*gc_x],
const int nx_, const int ny_, const int nx_, const int ny_,
const int boundary_conditions_) { const int boundary_conditions_) {
//Handle reflective boundary conditions //Handle reflective boundary conditions
if (getBCNorth(boundary_conditions_) == Reflective) { if (getBCNorth(boundary_conditions_) == Reflective) {
bcNorthReflective<block_width, block_height, ghost_cells, sign_y>(Q, nx_, ny_); bcNorthReflective<w, h, gc_x, gc_y, sign_y>(Q, nx_, ny_);
__syncthreads(); __syncthreads();
} }
if (getBCSouth(boundary_conditions_) == Reflective) { if (getBCSouth(boundary_conditions_) == Reflective) {
bcSouthReflective<block_width, block_height, ghost_cells, sign_y>(Q, nx_, ny_); bcSouthReflective<w, h, gc_x, gc_y, sign_y>(Q, nx_, ny_);
__syncthreads(); __syncthreads();
} }
if (getBCEast(boundary_conditions_) == Reflective) { if (getBCEast(boundary_conditions_) == Reflective) {
bcEastReflective<block_width, block_height, ghost_cells, sign_x>(Q, nx_, ny_); bcEastReflective<w, h, gc_x, gc_y, sign_x>(Q, nx_, ny_);
__syncthreads(); __syncthreads();
} }
if (getBCWest(boundary_conditions_) == Reflective) { if (getBCWest(boundary_conditions_) == Reflective) {
bcWestReflective<block_width, block_height, ghost_cells, sign_x>(Q, nx_, ny_); bcWestReflective<w, h, gc_x, gc_y, sign_x>(Q, nx_, ny_);
__syncthreads(); __syncthreads();
} }
} }
@ -204,9 +210,9 @@ inline __device__ int handleReflectiveBoundary(
/** /**
* Reads a block of data with ghost cells * Reads a block of data with ghost cells
*/ */
template<int block_width, int block_height, int ghost_cells, int sign_x, int sign_y> template<int w, int h, int gc_x, int gc_y, int sign_x, int sign_y>
inline __device__ void readBlock(float* ptr_, int pitch_, inline __device__ void readBlock(float* ptr_, int pitch_,
float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], float Q[h+2*gc_y][w+2*gc_x],
const int nx_, const int ny_, const int nx_, const int ny_,
const int boundary_conditions_) { const int boundary_conditions_) {
//Index of block within domain //Index of block within domain
@ -215,16 +221,16 @@ inline __device__ void readBlock(float* ptr_, int pitch_,
//Read into shared memory //Read into shared memory
//Loop over all variables //Loop over all variables
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) { for (int j=threadIdx.y; j<h+2*gc_y; j+=h) {
//Handle periodic boundary conditions here //Handle periodic boundary conditions here
int l = handlePeriodicBoundaryY<ghost_cells>(by + j, ny_, boundary_conditions_); int l = handlePeriodicBoundaryY<gc_y>(by + j, ny_, boundary_conditions_);
l = min(l, ny_+2*ghost_cells-1); l = min(l, ny_+2*gc_y-1);
float* row = (float*) ((char*) ptr_ + pitch_*l); float* row = (float*) ((char*) ptr_ + pitch_*l);
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) { for (int i=threadIdx.x; i<w+2*gc_x; i+=w) {
//Handle periodic boundary conditions here //Handle periodic boundary conditions here
int k = handlePeriodicBoundaryX<ghost_cells>(bx + i, nx_, boundary_conditions_); int k = handlePeriodicBoundaryX<gc_x>(bx + i, nx_, boundary_conditions_);
k = min(k, nx_+2*ghost_cells-1); k = min(k, nx_+2*gc_x-1);
//Read from global memory //Read from global memory
Q[j][i] = row[k]; Q[j][i] = row[k];
@ -232,7 +238,7 @@ inline __device__ void readBlock(float* ptr_, int pitch_,
} }
__syncthreads(); __syncthreads();
handleReflectiveBoundary<block_width, block_height, ghost_cells, sign_x, sign_y>(Q, nx_, ny_, boundary_conditions_); handleReflectiveBoundary<w, h, gc_x, gc_y, sign_x, sign_y>(Q, nx_, ny_, boundary_conditions_);
} }
@ -241,45 +247,68 @@ inline __device__ void readBlock(float* ptr_, int pitch_,
/** /**
* Writes a block of data to global memory for the shallow water equations. * Writes a block of data to global memory for the shallow water equations.
*/ */
template<int block_width, int block_height, int ghost_cells> template<int w, int h, int gc_x, int gc_y>
inline __device__ void writeBlock(float* ptr_, int pitch_, inline __device__ void writeBlock(float* ptr_, int pitch_,
float shmem[block_height+2*ghost_cells][block_width+2*ghost_cells], float shmem[h+2*gc_y][w+2*gc_x],
const int width, const int height, const int nx_, const int ny_,
int rk_step_, int rk_order_) { int rk_step_, int rk_order_) {
//Index of cell within domain //Index of cell within domain
const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells; const int ti = blockDim.x*blockIdx.x + threadIdx.x + gc_x;
const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells; const int tj = blockDim.y*blockIdx.y + threadIdx.y + gc_y;
//Only write internal cells //Only write internal cells
if (ti < width+ghost_cells && tj < height+ghost_cells) { if (ti < nx_+gc_x && tj < ny_+gc_y) {
//Index of thread within block //Index of thread within block
const int tx = threadIdx.x + ghost_cells; const int tx = threadIdx.x + gc_x;
const int ty = threadIdx.y + ghost_cells; const int ty = threadIdx.y + gc_y;
float* const row = (float*) ((char*) ptr_ + pitch_*tj); float* const row = (float*) ((char*) ptr_ + pitch_*tj);
//Handle runge-kutta timestepping here //Handle runge-kutta timestepping here
row[ti] = shmem[ty][tx]; row[ti] = shmem[ty][tx];
/**
* SSPRK1 (forward Euler)
* u^1 = u^n + dt*f(u^n)
*/
if (rk_order_ == 1) {
row[ti] = shmem[ty][tx];
}
/** /**
* SSPRK2 * SSPRK2
* u^1 = u^n + dt*f(u^n) * u^1 = u^n + dt*f(u^n)
* u^n+1 = 1/2*u^n + 1/2*(u^1 + dt*f(u^1)) * u^n+1 = 1/2*u^n + 1/2*(u^1 + dt*f(u^1))
* */
else if (rk_order_ == 2) {
if (rk_step_ == 0) {
row[ti] = shmem[ty][tx];
}
else if (rk_step_ == 1) {
row[ti] = 0.5f*row[ti] + 0.5f*shmem[ty][tx];
}
}
/**
* SSPRK3 * SSPRK3
* u^1 = u^n + dt*f(u^n) * u^1 = u^n + dt*f(u^n)
* u^2 = 3/4 * u^n + 1/4 * (u^1 + dt*f(u^1)) * u^2 = 3/4 * u^n + 1/4 * (u^1 + dt*f(u^1))
* u^n+1 = 1/3 * u^n + 2/3 * (u^2 + dt*f(u^2)) * u^n+1 = 1/3 * u^n + 2/3 * (u^2 + dt*f(u^2))
* FIXME: This is not correct now, need a temporary to hold intermediate step u^2
*/ */
else if (rk_order_ == 3) {
/* if (rk_step_ == 0) {
if (rk_order_ == 2 && rk_step_ == 1) { row[ti] = shmem[ty][tx];
row[ti] = 0.5f*(row[ti] + shmem[ty][tx]); }
else if (rk_step_ == 1) {
row[ti] = 0.75f*row[ti] + 0.25f*shmem[ty][tx];
}
else if (rk_step_ == 2) {
const float t = 1.0f / 3.0f; //Not representable in base 2
row[ti] = t*row[ti] + (1.0f-t)*shmem[ty][tx];
}
} }
else {
row[ti] = shmem[ty][tx];
}*/
} }
} }
@ -297,25 +326,26 @@ inline __device__ void writeBlock(float* ptr_, int pitch_,
// West boundary // West boundary
template<int block_width, int block_height, int ghost_cells, int sign> template<int w, int h, int gc_x, int gc_y, int sign>
__device__ void bcWestReflective(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) { __device__ void bcWestReflective(float Q[h+2*gc_y][w+2*gc_x],
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+= block_height) { const int nx_, const int ny_) {
const int i = threadIdx.x + ghost_cells; for (int j=threadIdx.y; j<h+2*gc_y; j+=h) {
const int i = threadIdx.x + gc_x;
const int ti = blockDim.x*blockIdx.x + i; const int ti = blockDim.x*blockIdx.x + i;
if (ti == ghost_cells) { if (gc_x >= 1 && ti == gc_x) {
Q[j][i-1] = sign*Q[j][i]; Q[j][i-1] = sign*Q[j][i];
} }
if (ghost_cells >= 2 && ti == ghost_cells + 1) { if (gc_x >= 2 && ti == gc_x + 1) {
Q[j][i-3] = sign*Q[j][i]; Q[j][i-3] = sign*Q[j][i];
} }
if (ghost_cells >= 3 && ti == ghost_cells + 2) { if (gc_x >= 3 && ti == gc_x + 2) {
Q[j][i-5] = sign*Q[j][i]; Q[j][i-5] = sign*Q[j][i];
} }
if (ghost_cells >= 4 && ti == ghost_cells + 3) { if (gc_x >= 4 && ti == gc_x + 3) {
Q[j][i-7] = sign*Q[j][i]; Q[j][i-7] = sign*Q[j][i];
} }
if (ghost_cells >= 5 && ti == ghost_cells + 4) { if (gc_x >= 5 && ti == gc_x + 4) {
Q[j][i-9] = sign*Q[j][i]; Q[j][i-9] = sign*Q[j][i];
} }
} }
@ -323,25 +353,26 @@ __device__ void bcWestReflective(float Q[block_height+2*ghost_cells][block_width
// East boundary // East boundary
template<int block_width, int block_height, int ghost_cells, int sign> template<int w, int h, int gc_x, int gc_y, int sign>
__device__ void bcEastReflective(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) { __device__ void bcEastReflective(float Q[h+2*gc_y][w+2*gc_x],
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+= block_height) { const int nx_, const int ny_) {
const int i = threadIdx.x + ghost_cells; for (int j=threadIdx.y; j<h+2*gc_y; j+=h) {
const int i = threadIdx.x + gc_x;
const int ti = blockDim.x*blockIdx.x + i; const int ti = blockDim.x*blockIdx.x + i;
if (ti == nx_ + ghost_cells - 1) { if (gc_x >= 1 && ti == nx_ + gc_x - 1) {
Q[j][i+1] = sign*Q[j][i]; Q[j][i+1] = sign*Q[j][i];
} }
if (ghost_cells >= 2 && ti == nx_ + ghost_cells - 2) { if (gc_x >= 2 && ti == nx_ + gc_x - 2) {
Q[j][i+3] = sign*Q[j][i]; Q[j][i+3] = sign*Q[j][i];
} }
if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 3) { if (gc_x >= 3 && ti == nx_ + gc_x - 3) {
Q[j][i+5] = sign*Q[j][i]; Q[j][i+5] = sign*Q[j][i];
} }
if (ghost_cells >= 4 && ti == nx_ + ghost_cells - 4) { if (gc_x >= 4 && ti == nx_ + gc_x - 4) {
Q[j][i+7] = sign*Q[j][i]; Q[j][i+7] = sign*Q[j][i];
} }
if (ghost_cells >= 5 && ti == nx_ + ghost_cells - 5) { if (gc_x >= 5 && ti == nx_ + gc_x - 5) {
Q[j][i+9] = sign*Q[j][i]; Q[j][i+9] = sign*Q[j][i];
} }
} }
@ -349,25 +380,26 @@ __device__ void bcEastReflective(float Q[block_height+2*ghost_cells][block_width
// South boundary // South boundary
template<int block_width, int block_height, int ghost_cells, int sign> template<int w, int h, int gc_x, int gc_y, int sign>
__device__ void bcSouthReflective(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) { __device__ void bcSouthReflective(float Q[h+2*gc_y][w+2*gc_x],
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+= block_width) { const int nx_, const int ny_) {
const int j = threadIdx.y + ghost_cells; for (int i=threadIdx.x; i<w+2*gc_x; i+=w) {
const int j = threadIdx.y + gc_y;
const int tj = blockDim.y*blockIdx.y + j; const int tj = blockDim.y*blockIdx.y + j;
if (tj == ghost_cells) { if (gc_y >= 1 && tj == gc_y) {
Q[j-1][i] = sign*Q[j][i]; Q[j-1][i] = sign*Q[j][i];
} }
if (ghost_cells >= 2 && tj == ghost_cells + 1) { if (gc_y >= 2 && tj == gc_y + 1) {
Q[j-3][i] = sign*Q[j][i]; Q[j-3][i] = sign*Q[j][i];
} }
if (ghost_cells >= 3 && tj == ghost_cells + 2) { if (gc_y >= 3 && tj == gc_y + 2) {
Q[j-5][i] = sign*Q[j][i]; Q[j-5][i] = sign*Q[j][i];
} }
if (ghost_cells >= 4 && tj == ghost_cells + 3) { if (gc_y >= 4 && tj == gc_y + 3) {
Q[j-7][i] = sign*Q[j][i]; Q[j-7][i] = sign*Q[j][i];
} }
if (ghost_cells >= 5 && tj == ghost_cells + 4) { if (gc_y >= 5 && tj == gc_y + 4) {
Q[j-9][i] = sign*Q[j][i]; Q[j-9][i] = sign*Q[j][i];
} }
} }
@ -377,25 +409,25 @@ __device__ void bcSouthReflective(float Q[block_height+2*ghost_cells][block_widt
// North boundary // North boundary
template<int block_width, int block_height, int ghost_cells, int sign> template<int w, int h, int gc_x, int gc_y, int sign>
__device__ void bcNorthReflective(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) { __device__ void bcNorthReflective(float Q[h+2*gc_y][w+2*gc_x], const int nx_, const int ny_) {
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+= block_width) { for (int i=threadIdx.x; i<w+2*gc_x; i+=w) {
const int j = threadIdx.y + ghost_cells; const int j = threadIdx.y + gc_y;
const int tj = blockDim.y*blockIdx.y + j; const int tj = blockDim.y*blockIdx.y + j;
if (tj == ny_ + ghost_cells - 1) { if (gc_y >= 1 && tj == ny_ + gc_y - 1) {
Q[j+1][i] = sign*Q[j][i]; Q[j+1][i] = sign*Q[j][i];
} }
if (ghost_cells >= 2 && tj == ny_ + ghost_cells - 2) { if (gc_y >= 2 && tj == ny_ + gc_y - 2) {
Q[j+3][i] = sign*Q[j][i]; Q[j+3][i] = sign*Q[j][i];
} }
if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 3) { if (gc_y >= 3 && tj == ny_ + gc_y - 3) {
Q[j+5][i] = sign*Q[j][i]; Q[j+5][i] = sign*Q[j][i];
} }
if (ghost_cells >= 4 && tj == ny_ + ghost_cells - 4) { if (gc_y >= 4 && tj == ny_ + gc_y - 4) {
Q[j+7][i] = sign*Q[j][i]; Q[j+7][i] = sign*Q[j][i];
} }
if (ghost_cells >= 5 && tj == ny_ + ghost_cells - 5) { if (gc_y >= 5 && tj == ny_ + gc_y - 5) {
Q[j+9][i] = sign*Q[j][i]; Q[j+9][i] = sign*Q[j][i];
} }
} }
@ -422,13 +454,13 @@ __device__ void bcNorthReflective(float Q[block_height+2*ghost_cells][block_widt
template<int block_width, int block_height, int ghost_cells, int vars> template<int w, int h, int gc_x, int gc_y, int vars>
__device__ void evolveF(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], __device__ void evolveF(float Q[vars][h+2*gc_y][w+2*gc_x],
float F[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], float F[vars][h+2*gc_y][w+2*gc_x],
const float dx_, const float dt_) { const float dx_, const float dt_) {
for (int var=0; var < vars; ++var) { for (int var=0; var < vars; ++var) {
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) { for (int j=threadIdx.y; j<h+2*gc_y; j+=h) {
for (int i=threadIdx.x+ghost_cells; i<block_width+ghost_cells; i+=block_width) { for (int i=threadIdx.x+gc_x; i<w+gc_x; i+=w) {
Q[var][j][i] = Q[var][j][i] + (F[var][j][i-1] - F[var][j][i]) * dt_ / dx_; Q[var][j][i] = Q[var][j][i] + (F[var][j][i-1] - F[var][j][i]) * dt_ / dx_;
} }
} }
@ -443,13 +475,13 @@ __device__ void evolveF(float Q[vars][block_height+2*ghost_cells][block_width+2*
/** /**
* Evolves the solution in time along the y axis (dimensional splitting) * Evolves the solution in time along the y axis (dimensional splitting)
*/ */
template<int block_width, int block_height, int ghost_cells, int vars> template<int w, int h, int gc_x, int gc_y, int vars>
__device__ void evolveG(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], __device__ void evolveG(float Q[vars][h+2*gc_y][w+2*gc_x],
float G[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], float G[vars][h+2*gc_y][w+2*gc_x],
const float dy_, const float dt_) { const float dy_, const float dt_) {
for (int var=0; var < vars; ++var) { for (int var=0; var < vars; ++var) {
for (int j=threadIdx.y+ghost_cells; j<block_height+ghost_cells; j+=block_height) { for (int j=threadIdx.y+gc_y; j<h+gc_y; j+=h) {
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) { for (int i=threadIdx.x; i<w+2*gc_x; i+=w) {
Q[var][j][i] = Q[var][j][i] + (G[var][j-1][i] - G[var][j][i]) * dt_ / dy_; Q[var][j][i] = Q[var][j][i] + (G[var][j-1][i] - G[var][j][i]) * dt_ / dy_;
} }
} }
@ -478,6 +510,55 @@ __device__ void memset(float Q[vars][shmem_height][shmem_width], float value) {
template <unsigned int threads>
__device__ void reduce_max(float* data, unsigned int n) {
__shared__ float sdata[threads];
unsigned int tid = threadIdx.x;
//Reduce to "threads" elements
sdata[tid] = FLT_MIN;
for (unsigned int i=tid; i<n; i += threads) {
sdata[tid] = max(sdata[tid], dt_ctx.L[i]);
}
__syncthreads();
//Now, reduce all elements into a single element
if (threads >= 512) {
if (tid < 256) {
sdata[tid] = max(sdata[tid], sdata[tid + 256]);
}
__syncthreads();
}
if (threads >= 256) {
if (tid < 128) {
sdata[tid] = max(sdata[tid], sdata[tid + 128]);
}
__syncthreads();
}
if (threads >= 128) {
if (tid < 64) {
sdata[tid] = max(sdata[tid], sdata[tid + 64]);
}
__syncthreads();
}
if (tid < 32) {
volatile float* sdata_volatile = sdata;
if (threads >= 64) {
sdata_volatile[tid] = max(sdata_volatile[tid], sdata_volatile[tid + 32]);
}
if (tid < 16) {
if (threads >= 32) sdata_volatile[tid] = max(sdata_volatile[tid], sdata_volatile[tid + 16]);
if (threads >= 16) sdata_volatile[tid] = max(sdata_volatile[tid], sdata_volatile[tid + 8]);
if (threads >= 8) sdata_volatile[tid] = max(sdata_volatile[tid], sdata_volatile[tid + 4]);
if (threads >= 4) sdata_volatile[tid] = max(sdata_volatile[tid], sdata_volatile[tid + 2]);
if (threads >= 2) sdata_volatile[tid] = max(sdata_volatile[tid], sdata_volatile[tid + 1]);
}
if (tid == 0) {
return sdata_volatile[0];
}
}
}

20
GPUSimulators/cuda/limiters.h
View File

@ -46,14 +46,14 @@ __device__ __inline__ float minmodSlope(float left, float center, float right, f
/** /**
* Reconstructs a minmod slope for a whole block along the abscissa * Reconstructs a minmod slope for a whole block along the abscissa
*/ */
template<int block_width, int block_height, int ghost_cells, int vars> template<int w, int h, int gc_x, int gc_y, int vars>
__device__ void minmodSlopeX(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], __device__ void minmodSlopeX(float Q[vars][h+2*gc_y][w+2*gc_x],
float Qx[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], float Qx[vars][h+2*gc_y][w+2*gc_x],
const float theta_) { const float theta_) {
//Reconstruct slopes along x axis //Reconstruct slopes along x axis
for (int p=0; p<vars; ++p) { for (int p=0; p<vars; ++p) {
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) { for (int j=threadIdx.y; j<h+2*gc_y; j+=h) {
for (int i=threadIdx.x+1; i<block_width+3; i+=block_width) { for (int i=threadIdx.x+1; i<w+2*gc_x-1; i+=w) {
Qx[p][j][i] = minmodSlope(Q[p][j][i-1], Q[p][j][i], Q[p][j][i+1], theta_); Qx[p][j][i] = minmodSlope(Q[p][j][i-1], Q[p][j][i], Q[p][j][i+1], theta_);
} }
} }
@ -64,14 +64,14 @@ __device__ void minmodSlopeX(float Q[vars][block_height+2*ghost_cells][block_wid
/** /**
* Reconstructs a minmod slope for a whole block along the ordinate * Reconstructs a minmod slope for a whole block along the ordinate
*/ */
template<int block_width, int block_height, int ghost_cells, int vars> template<int w, int h, int gc_x, int gc_y, int vars>
__device__ void minmodSlopeY(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], __device__ void minmodSlopeY(float Q[vars][h+2*gc_y][w+2*gc_x],
float Qy[vars][block_height+2*ghost_cells][block_width+2*ghost_cells], float Qy[vars][h+2*gc_y][w+2*gc_x],
const float theta_) { const float theta_) {
//Reconstruct slopes along y axis //Reconstruct slopes along y axis
for (int p=0; p<vars; ++p) { for (int p=0; p<vars; ++p) {
for (int j=threadIdx.y+1; j<block_height+3; j+=block_height) { for (int j=threadIdx.y+1; j<h+2*gc_y-1; j+=h) {
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) { for (int i=threadIdx.x; i<w+2*gc_x; i+=w) {
Qy[p][j][i] = minmodSlope(Q[p][j-1][i], Q[p][j][i], Q[p][j+1][i], theta_); Qy[p][j][i] = minmodSlope(Q[p][j-1][i], Q[p][j][i], Q[p][j+1][i], theta_);
} }
} }