From 2d8858e7e6cdd85318266c3d89cac311e8c78a3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20R=2E=20Brodtkorb?= Date: Wed, 31 Oct 2018 15:34:54 +0100 Subject: [PATCH] Refactoring --- .../{HLL2Euler.py => EE2D_KP07_dimsplit.py} | 43 +-- GPUSimulators/FORCE.py | 5 +- GPUSimulators/HLL.py | 5 +- GPUSimulators/HLL2.py | 6 +- GPUSimulators/KP07.py | 7 +- GPUSimulators/KP07_dimsplit.py | 9 +- GPUSimulators/LxF.py | 5 +- GPUSimulators/Simulator.py | 3 - GPUSimulators/WAF.py | 6 +- GPUSimulators/cuda/EE2D_KP07_dimsplit.cu | 244 ++++++++++++++++++ GPUSimulators/cuda/EulerCommon.h | 8 +- .../cuda/{SWE_FORCE.cu => SWE2D_FORCE.cu} | 10 +- .../cuda/{SWE_HLL.cu => SWE2D_HLL.cu} | 10 +- .../cuda/{SWE_HLL2.cu => SWE2D_HLL2.cu} | 20 +- .../cuda/{SWE_KP07.cu => SWE2D_KP07.cu} | 0 ...P07_dimsplit.cu => SWE2D_KP07_dimsplit.cu} | 27 +- .../cuda/{SWE_LxF.cu => SWE2D_LxF.cu} | 0 .../cuda/{SWE_WAF.cu => SWE2D_WAF.cu} | 20 +- GPUSimulators/cuda/common.h | 132 +++------- 19 files changed, 402 insertions(+), 158 deletions(-) rename GPUSimulators/{HLL2Euler.py => EE2D_KP07_dimsplit.py} (76%) create mode 100644 GPUSimulators/cuda/EE2D_KP07_dimsplit.cu rename GPUSimulators/cuda/{SWE_FORCE.cu => SWE2D_FORCE.cu} (94%) rename GPUSimulators/cuda/{SWE_HLL.cu => SWE2D_HLL.cu} (94%) rename GPUSimulators/cuda/{SWE_HLL2.cu => SWE2D_HLL2.cu} (93%) rename GPUSimulators/cuda/{SWE_KP07.cu => SWE2D_KP07.cu} (100%) rename GPUSimulators/cuda/{SWE_KP07_dimsplit.cu => SWE2D_KP07_dimsplit.cu} (91%) rename GPUSimulators/cuda/{SWE_LxF.cu => SWE2D_LxF.cu} (100%) rename GPUSimulators/cuda/{SWE_WAF.cu => SWE2D_WAF.cu} (90%) diff --git a/GPUSimulators/HLL2Euler.py b/GPUSimulators/EE2D_KP07_dimsplit.py similarity index 76% rename from GPUSimulators/HLL2Euler.py rename to GPUSimulators/EE2D_KP07_dimsplit.py index 9346d9c..c98136f 100644 --- a/GPUSimulators/HLL2Euler.py +++ b/GPUSimulators/EE2D_KP07_dimsplit.py @@ -20,8 +20,8 @@ along with this program. If not, see . """ #Import packages we need -import numpy as np from GPUSimulators import Simulator, Common +import numpy as np @@ -34,7 +34,7 @@ from GPUSimulators import Simulator, Common """ Class that solves the SW equations using the Forward-Backward linear scheme """ -class HLL2Euler (Simulator.BaseSimulator): +class EE2D_KP07_dimsplit (Simulator.BaseSimulator): """ Initialization routine @@ -47,31 +47,40 @@ class HLL2Euler (Simulator.BaseSimulator): dx: Grid cell spacing along x-axis dy: Grid cell spacing along y-axis dt: Size of each timestep - g: Gravitational accelleration + gamma: Gas constant + p: pressure """ def __init__(self, \ context, \ rho, rho_u, rho_v, E, \ nx, ny, \ dx, dy, dt, \ - g, \ - theta=1.8, \ + gamma, \ + theta=1.3, \ block_width=16, block_height=16): # Call super constructor super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ - block_width, block_height); + block_width, block_height) + self.gamma = np.float32(gamma) self.theta = np.float32(theta) #Get kernels - self.kernel = context.get_prepared_kernel("HLL2_kernel.cu", "HLL2Kernel", \ - "iifffffiPiPiPiPiPiPi", \ - BLOCK_WIDTH=self.local_size[0], \ - BLOCK_HEIGHT=self.local_size[1]) + #Get kernels + self.kernel = context.get_prepared_kernel("cuda/EE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \ + "iifffffiPiPiPiPiPiPiPiPi", \ + defines={ + 'BLOCK_WIDTH': self.block_size[0], + 'BLOCK_HEIGHT': self.block_size[1] + }, \ + compile_args={ + 'no_extern_c': True, + 'options': ["--use_fast_math"], + }, \ + jit_compile_args={}) #Create data by uploading to device self.u0 = Common.ArakawaA2D(self.stream, \ @@ -90,10 +99,10 @@ class HLL2Euler (Simulator.BaseSimulator): return self.stepDimsplitXY(dt) def stepDimsplitXY(self, dt): - self.kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \ + self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.nx, self.ny, \ self.dx, self.dy, dt, \ - self.g, \ + self.gamma, \ self.theta, \ np.int32(0), \ self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ @@ -102,16 +111,16 @@ class HLL2Euler (Simulator.BaseSimulator): self.u0[3].data.gpudata, self.u0[3].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[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.u0, self.u1 = self.u1, self.u0 self.t += dt def stepDimsplitYX(self, dt): - self.kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \ + self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \ self.nx, self.ny, \ self.dx, self.dy, dt, \ - self.g, \ + self.gamma, \ self.theta, \ np.int32(1), \ self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ @@ -120,7 +129,7 @@ class HLL2Euler (Simulator.BaseSimulator): self.u0[3].data.gpudata, self.u0[3].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[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.u0, self.u1 = self.u1, self.u0 self.t += dt diff --git a/GPUSimulators/FORCE.py b/GPUSimulators/FORCE.py index e88224c..bd8074c 100644 --- a/GPUSimulators/FORCE.py +++ b/GPUSimulators/FORCE.py @@ -22,6 +22,7 @@ along with this program. If not, see . #Import packages we need from GPUSimulators import Simulator, Common +import numpy as np @@ -62,11 +63,11 @@ class FORCE (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height); + self.g = np.float32(g) #Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_FORCE.cu", "FORCEKernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_FORCE.cu", "FORCEKernel", \ "iiffffPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/HLL.py b/GPUSimulators/HLL.py index bc77ff8..11e5219 100644 --- a/GPUSimulators/HLL.py +++ b/GPUSimulators/HLL.py @@ -21,6 +21,7 @@ along with this program. If not, see . #Import packages we need from GPUSimulators import Simulator, Common +import numpy as np @@ -57,11 +58,11 @@ class HLL (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height); + self.g = np.float32(g) #Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_HLL.cu", "HLLKernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_HLL.cu", "HLLKernel", \ "iiffffPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/HLL2.py b/GPUSimulators/HLL2.py index f4e141b..e45291d 100644 --- a/GPUSimulators/HLL2.py +++ b/GPUSimulators/HLL2.py @@ -20,8 +20,8 @@ along with this program. If not, see . """ #Import packages we need -import numpy as np from GPUSimulators import Simulator, Common +import numpy as np @@ -61,13 +61,13 @@ class HLL2 (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height); + self.g = np.float32(g) self.theta = np.float32(theta) #Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_HLL2.cu", "HLL2Kernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_HLL2.cu", "HLL2Kernel", \ "iifffffiPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/KP07.py b/GPUSimulators/KP07.py index 46ad8f4..75e4a11 100644 --- a/GPUSimulators/KP07.py +++ b/GPUSimulators/KP07.py @@ -25,8 +25,8 @@ along with this program. If not, see . """ #Import packages we need -import numpy as np from GPUSimulators import Simulator, Common +import numpy as np @@ -62,13 +62,12 @@ class KP07 (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height); - + self.g = np.float32(g) self.theta = np.float32(theta) #Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_KP07.cu", "KP07Kernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_KP07.cu", "KP07Kernel", \ "iifffffiPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/KP07_dimsplit.py b/GPUSimulators/KP07_dimsplit.py index 92e7fa5..4fd080e 100644 --- a/GPUSimulators/KP07_dimsplit.py +++ b/GPUSimulators/KP07_dimsplit.py @@ -25,8 +25,8 @@ along with this program. If not, see . """ #Import packages we need -import numpy as np from GPUSimulators import Simulator, Common +import numpy as np @@ -62,13 +62,12 @@ class KP07_dimsplit (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ - block_width, block_height); - + block_width, block_height) + self.g = np.float32(g) self.theta = np.float32(theta) #Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_KP07_dimsplit.cu", "KP07DimsplitKernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \ "iifffffiPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/LxF.py b/GPUSimulators/LxF.py index 33ab080..c8aa554 100644 --- a/GPUSimulators/LxF.py +++ b/GPUSimulators/LxF.py @@ -22,6 +22,7 @@ along with this program. If not, see . #Import packages we need from GPUSimulators import Simulator, Common, CudaContext +import numpy as np @@ -58,11 +59,11 @@ class LxF (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height); + self.g = np.float32(g) # Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_LxF.cu", "LxFKernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_LxF.cu", "LxFKernel", \ "iiffffPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/Simulator.py b/GPUSimulators/Simulator.py index 28fcd7a..3bf7abf 100644 --- a/GPUSimulators/Simulator.py +++ b/GPUSimulators/Simulator.py @@ -44,13 +44,11 @@ class BaseSimulator: dx: Grid cell spacing along x-axis (20 000 m) dy: Grid cell spacing along y-axis (20 000 m) dt: Size of each timestep (90 s) - g: Gravitational accelleration (9.81 m/s^2) """ def __init__(self, \ context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height): #Get logger self.logger = logging.getLogger(__name__ + "." + self.__class__.__name__) @@ -64,7 +62,6 @@ class BaseSimulator: self.dx = np.float32(dx) self.dy = np.float32(dy) self.dt = np.float32(dt) - self.g = np.float32(g) #Handle autotuning block size if (self.context.autotuner): diff --git a/GPUSimulators/WAF.py b/GPUSimulators/WAF.py index c22c7f6..d9ae132 100644 --- a/GPUSimulators/WAF.py +++ b/GPUSimulators/WAF.py @@ -21,8 +21,8 @@ along with this program. If not, see . """ #Import packages we need -import numpy as np from GPUSimulators import Simulator, Common +import numpy as np @@ -57,11 +57,11 @@ class WAF (Simulator.BaseSimulator): super().__init__(context, \ nx, ny, \ dx, dy, dt, \ - g, \ block_width, block_height); + self.g = np.float32(g) #Get kernels - self.kernel = context.get_prepared_kernel("cuda/SWE_WAF.cu", "WAFKernel", \ + self.kernel = context.get_prepared_kernel("cuda/SWE2D_WAF.cu", "WAFKernel", \ "iiffffiPiPiPiPiPiPi", \ defines={ 'BLOCK_WIDTH': self.block_size[0], diff --git a/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu new file mode 100644 index 0000000..c404ed2 --- /dev/null +++ b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu @@ -0,0 +1,244 @@ + /* +This kernel implements the Central Upwind flux function to +solve the Euler equations + +Copyright (C) 2018 SINTEF Digital + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +*/ + + + +#include "common.h" +#include "EulerCommon.h" +#include "limiters.h" + + +__device__ +void computeFluxF(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], + float Qx[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2], + float F[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1], + const float gamma_, const float dx_, const float dt_) { + //Index of thread within block + const int tx = threadIdx.x; + const int ty = threadIdx.y; + + { + int j=ty; + const int l = j + 2; //Skip ghost cells + for (int i=tx; i( rho0_ptr_, rho0_pitch_, Q[0], nx_+2, ny_+2); + readBlock(rho_u0_ptr_, rho_u0_pitch_, Q[1], nx_+2, ny_+2); + readBlock(rho_v0_ptr_, rho_v0_pitch_, Q[2], nx_+2, ny_+2); + readBlock( E0_ptr_, E0_pitch_, Q[3], nx_+2, ny_+2); + __syncthreads(); + + + //Fix boundary conditions + noFlowBoundary(Q[0], nx_, ny_); + noFlowBoundary(Q[1], nx_, ny_); + noFlowBoundary(Q[2], nx_, ny_); + noFlowBoundary(Q[3], nx_, ny_); + __syncthreads(); + + + //Step 0 => evolve x first, then y + if (step_ == 0) { + //Compute fluxes along the x axis and evolve + minmodSlopeX(Q, Qx, theta_); + __syncthreads(); + computeFluxF(Q, Qx, F, gamma_, dx_, dt_); + __syncthreads(); + evolveF2(Q, F, nx_, ny_, dx_, dt_); + __syncthreads(); + + //Set boundary conditions + noFlowBoundary(Q[0], nx_, ny_); + noFlowBoundary(Q[1], nx_, ny_); + noFlowBoundary(Q[2], nx_, ny_); + noFlowBoundary(Q[3], nx_, ny_); + __syncthreads(); + + //Compute fluxes along the y axis and evolve + minmodSlopeY(Q, Qx, theta_); + __syncthreads(); + computeFluxG(Q, Qx, F, gamma_, dy_, dt_); + __syncthreads(); + evolveG2(Q, F, nx_, ny_, dy_, dt_); + __syncthreads(); + } + //Step 1 => evolve y first, then x + else { + //Compute fluxes along the y axis and evolve + minmodSlopeY(Q, Qx, theta_); + __syncthreads(); + computeFluxG(Q, Qx, F, gamma_, dy_, dt_); + __syncthreads(); + evolveG2(Q, F, nx_, ny_, dy_, dt_); + __syncthreads(); + + //Set boundary conditions + noFlowBoundary(Q[0], nx_, ny_); + noFlowBoundary(Q[1], nx_, ny_); + noFlowBoundary(Q[2], nx_, ny_); + noFlowBoundary(Q[3], nx_, ny_); + __syncthreads(); + + //Compute fluxes along the x axis and evolve + minmodSlopeX(Q, Qx, theta_); + __syncthreads(); + computeFluxF(Q, Qx, F, gamma_, dx_, dt_); + __syncthreads(); + evolveF2(Q, F, nx_, ny_, dx_, dt_); + __syncthreads(); + } + + + // Write to main memory for all internal cells + writeBlock( rho1_ptr_, rho1_pitch_, Q[0], nx_, ny_); + writeBlock(rho_u1_ptr_, rho_u1_pitch_, Q[1], nx_, ny_); + writeBlock(rho_v1_ptr_, rho_v1_pitch_, Q[2], nx_, ny_); + writeBlock( E1_ptr_, E1_pitch_, Q[3], nx_, ny_); +} + +} // extern "C" \ No newline at end of file diff --git a/GPUSimulators/cuda/EulerCommon.h b/GPUSimulators/cuda/EulerCommon.h index 47b1eb2..a748e0f 100644 --- a/GPUSimulators/cuda/EulerCommon.h +++ b/GPUSimulators/cuda/EulerCommon.h @@ -66,15 +66,15 @@ __device__ float4 CentralUpwindFlux(const float4 Qm, float4 Qp, const float gamm const float Pp = pressure(Qp, gamma); const float4 Fp = F_func(Qp, Pp); const float up = Qp.y / Qp.x; // rho*u / rho - const float cp = sqrt(gamma*P*Qp.x); // sqrt(gamma*P/rho) + const float cp = sqrt(gamma*Pp*Qp.x); // sqrt(gamma*P/rho) const float Pm = pressure(Qm, gamma); - const float3 Fm = F_func(Qm, Pm); + const float4 Fm = F_func(Qm, Pm); const float um = Qm.y / Qm.x; // rho*u / rho - const float cm = sqrt(gamma*P/Qm.x); // sqrt(g*h) + const float cm = sqrt(gamma*Pm/Qm.x); // sqrt(g*h) const float am = min(min(um-cm, up-cp), 0.0f); // largest negative wave speed const float ap = max(max(um+cm, up+cp), 0.0f); // largest positive wave speed - return ((ap*Fm - am*Fp) + ap*am*(Qp-Qm))/(ap-am); + return ((ap*Fm - am*Fp) + ap*am*(Qp-Qm))/(ap-am); } \ No newline at end of file diff --git a/GPUSimulators/cuda/SWE_FORCE.cu b/GPUSimulators/cuda/SWE2D_FORCE.cu similarity index 94% rename from GPUSimulators/cuda/SWE_FORCE.cu rename to GPUSimulators/cuda/SWE2D_FORCE.cu index adab937..0762029 100644 --- a/GPUSimulators/cuda/SWE_FORCE.cu +++ b/GPUSimulators/cuda/SWE2D_FORCE.cu @@ -136,7 +136,10 @@ __global__ void FORCEKernel( //Compute flux along x, and evolve computeFluxF(Q, F, g_, dx_, dt_); __syncthreads(); - evolveF1(Q, F, nx_, ny_, dx_, dt_); + + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); //Set boundary conditions @@ -148,7 +151,10 @@ __global__ void FORCEKernel( //Compute flux along y, and evolve computeFluxG(Q, F, g_, dy_, dt_); __syncthreads(); - evolveG1(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); //Write to main memory diff --git a/GPUSimulators/cuda/SWE_HLL.cu b/GPUSimulators/cuda/SWE2D_HLL.cu similarity index 94% rename from GPUSimulators/cuda/SWE_HLL.cu rename to GPUSimulators/cuda/SWE2D_HLL.cu index e34658c..8b0c2cb 100644 --- a/GPUSimulators/cuda/SWE_HLL.cu +++ b/GPUSimulators/cuda/SWE2D_HLL.cu @@ -144,7 +144,10 @@ __global__ void HLLKernel( //Compute F flux computeFluxF(Q, F, g_); __syncthreads(); - evolveF1(Q, F, nx_, ny_, dx_, dt_); + + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); //Set boundary conditions @@ -156,7 +159,10 @@ __global__ void HLLKernel( //Compute G flux computeFluxG(Q, F, g_); __syncthreads(); - evolveG1(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); // Write to main memory for all internal cells diff --git a/GPUSimulators/cuda/SWE_HLL2.cu b/GPUSimulators/cuda/SWE2D_HLL2.cu similarity index 93% rename from GPUSimulators/cuda/SWE_HLL2.cu rename to GPUSimulators/cuda/SWE2D_HLL2.cu index dc63d66..40910cc 100644 --- a/GPUSimulators/cuda/SWE_HLL2.cu +++ b/GPUSimulators/cuda/SWE2D_HLL2.cu @@ -184,7 +184,10 @@ __global__ void HLL2Kernel( __syncthreads(); computeFluxF(Q, Qx, F, g_, dx_, dt_); __syncthreads(); - evolveF2(Q, F, nx_, ny_, dx_, dt_); + + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); //Set boundary conditions @@ -198,7 +201,10 @@ __global__ void HLL2Kernel( __syncthreads(); computeFluxG(Q, Qx, F, g_, dy_, dt_); __syncthreads(); - evolveG2(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); } //Step 1 => evolve y first, then x @@ -208,7 +214,10 @@ __global__ void HLL2Kernel( __syncthreads(); computeFluxG(Q, Qx, F, g_, dy_, dt_); __syncthreads(); - evolveG2(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); //Set boundary conditions @@ -222,7 +231,10 @@ __global__ void HLL2Kernel( __syncthreads(); computeFluxF(Q, Qx, F, g_, dx_, dt_); __syncthreads(); - evolveF2(Q, F, nx_, ny_, dx_, dt_); + + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); } diff --git a/GPUSimulators/cuda/SWE_KP07.cu b/GPUSimulators/cuda/SWE2D_KP07.cu similarity index 100% rename from GPUSimulators/cuda/SWE_KP07.cu rename to GPUSimulators/cuda/SWE2D_KP07.cu diff --git a/GPUSimulators/cuda/SWE_KP07_dimsplit.cu b/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu similarity index 91% rename from GPUSimulators/cuda/SWE_KP07_dimsplit.cu rename to GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu index 74381d0..2957dcc 100644 --- a/GPUSimulators/cuda/SWE_KP07_dimsplit.cu +++ b/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu @@ -181,7 +181,9 @@ __global__ void KP07DimsplitKernel( __syncthreads(); computeFluxF(Q, Qx, F, g_, dx_, dt_); __syncthreads(); - evolveF2(Q, F, nx_, ny_, dx_, dt_); + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); //Set boundary conditions @@ -190,12 +192,18 @@ __global__ void KP07DimsplitKernel( noFlowBoundary(Q[2], nx_, ny_); __syncthreads(); + + //Compute fluxes along the y axis and evolve minmodSlopeY(Q, Qx, theta_); __syncthreads(); + computeFluxG(Q, Qx, F, g_, dy_, dt_); __syncthreads(); - evolveG2(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); } //Step 1 => evolve y first, then x @@ -205,7 +213,10 @@ __global__ void KP07DimsplitKernel( __syncthreads(); computeFluxG(Q, Qx, F, g_, dy_, dt_); __syncthreads(); - evolveG2(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); //Set boundary conditions @@ -219,15 +230,17 @@ __global__ void KP07DimsplitKernel( __syncthreads(); computeFluxF(Q, Qx, F, g_, dx_, dt_); __syncthreads(); - evolveF2(Q, F, nx_, ny_, dx_, dt_); + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); } // Write to main memory for all internal cells - writeBlock( h1_ptr_, h1_pitch_, Q[0], nx_, ny_); - writeBlock(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_); - writeBlock(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_); + writeBlock( h1_ptr_, h1_pitch_, Q[0], nx_, ny_); + writeBlock(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_); + writeBlock(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_); } } // extern "C" \ No newline at end of file diff --git a/GPUSimulators/cuda/SWE_LxF.cu b/GPUSimulators/cuda/SWE2D_LxF.cu similarity index 100% rename from GPUSimulators/cuda/SWE_LxF.cu rename to GPUSimulators/cuda/SWE2D_LxF.cu diff --git a/GPUSimulators/cuda/SWE_WAF.cu b/GPUSimulators/cuda/SWE2D_WAF.cu similarity index 90% rename from GPUSimulators/cuda/SWE_WAF.cu rename to GPUSimulators/cuda/SWE2D_WAF.cu index 62a82e2..1b80e4d 100644 --- a/GPUSimulators/cuda/SWE_WAF.cu +++ b/GPUSimulators/cuda/SWE2D_WAF.cu @@ -161,7 +161,10 @@ __global__ void WAFKernel( //Compute fluxes along the x axis and evolve computeFluxF(Q, F, g_, dx_, dt_); __syncthreads(); - evolveF2(Q, F, nx_, ny_, dx_, dt_); + + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); //Fix boundary conditions @@ -173,7 +176,10 @@ __global__ void WAFKernel( //Compute fluxes along the y axis and evolve computeFluxG(Q, F, g_, dy_, dt_); __syncthreads(); - evolveG2(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); } //Step 1 => evolve y first, then x @@ -181,7 +187,10 @@ __global__ void WAFKernel( //Compute fluxes along the y axis and evolve computeFluxG(Q, F, g_, dy_, dt_); __syncthreads(); - evolveG2(Q, F, nx_, ny_, dy_, dt_); + + evolveG(Q[0], F[0], dy_, dt_); + evolveG(Q[1], F[1], dy_, dt_); + evolveG(Q[2], F[2], dy_, dt_); __syncthreads(); //Fix boundary conditions @@ -193,7 +202,10 @@ __global__ void WAFKernel( //Compute fluxes along the x axis and evolve computeFluxF(Q, F, g_, dx_, dt_); __syncthreads(); - evolveF2(Q, F, nx_, ny_, dx_, dt_); + + evolveF(Q[0], F[0], dx_, dt_); + evolveF(Q[1], F[1], dx_, dt_); + evolveF(Q[2], F[2], dx_, dt_); __syncthreads(); } diff --git a/GPUSimulators/cuda/common.h b/GPUSimulators/cuda/common.h index faf961d..9fe08ff 100644 --- a/GPUSimulators/cuda/common.h +++ b/GPUSimulators/cuda/common.h @@ -44,6 +44,28 @@ inline __device__ float3 operator+(const float3 a, const float3 b) { return make_float3(a.x+b.x, a.y+b.y, a.z+b.z); } +/** + * Float4 operators + */ +inline __device__ float4 operator*(const float a, const float4 b) { + return make_float4(a*b.x, a*b.y, a*b.z, a*b.w); +} + +inline __device__ float4 operator/(const float4 a, const float b) { + return make_float4(a.x/b, a.y/b, a.z/b, a.w/b); +} + +inline __device__ float4 operator-(const float4 a, const float4 b) { + return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w); +} + +inline __device__ float4 operator+(const float4 a, const float4 b) { + return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w); +} + + + + inline __device__ __host__ float clamp(const float f, const float a, const float b) { return fmaxf(a, fminf(f, b)); } @@ -223,62 +245,23 @@ __device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2 - - - -/** - * Evolves the solution in time along the x axis (dimensional splitting) - */ -__device__ void evolveF1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2], - float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1], - const int nx_, const int ny_, +template +__device__ void evolveF(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], + float F[block_height+1][block_width+1], const float dx_, const float dt_) { //Index of thread within block const int tx = threadIdx.x; const int ty = threadIdx.y; + + const int i = tx + ghost_cells; //Skip local ghost cells + const int j = ty + ghost_cells; //Index of cell within domain - const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1 - const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1; + //const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells; //Skip global ghost cells, i.e., +1 + //const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells; + //if (ti > ghost_cells-1 && ti < nx_+ghost_cells && tj > ghost_cells-1 && tj < ny_+ghost_cells) { + Q[j][i] = Q[j][i] + (F[ty][tx] - F[ty][tx+1]) * dt_ / dx_; - if (ti > 0 && ti < nx_+1 && tj > 0 && tj < ny_+1) { - const int i = tx + 1; //Skip local ghost cells, i.e., +1 - const int j = ty + 1; - - Q[0][j][i] = Q[0][j][i] + (F[0][ty][tx] - F[0][ty][tx+1]) * dt_ / dx_; - Q[1][j][i] = Q[1][j][i] + (F[1][ty][tx] - F[1][ty][tx+1]) * dt_ / dx_; - Q[2][j][i] = Q[2][j][i] + (F[2][ty][tx] - F[2][ty][tx+1]) * dt_ / dx_; - } -} - - - - - - -/** - * Evolves the solution in time along the x axis (dimensional splitting) - */ -__device__ void evolveF2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], - float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1], - const int nx_, const int ny_, - const float dx_, const float dt_) { - //Index of thread within block - const int tx = threadIdx.x; - const int ty = threadIdx.y; - - //Index of cell within domain - const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2 - const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2; - - if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) { - const int i = tx + 2; //Skip local ghost cells, i.e., +1 - const int j = ty + 2; - - Q[0][j][i] = Q[0][j][i] + (F[0][ty][tx] - F[0][ty][tx+1]) * dt_ / dx_; - Q[1][j][i] = Q[1][j][i] + (F[1][ty][tx] - F[1][ty][tx+1]) * dt_ / dx_; - Q[2][j][i] = Q[2][j][i] + (F[2][ty][tx] - F[2][ty][tx+1]) * dt_ / dx_; - } } @@ -289,57 +272,18 @@ __device__ void evolveF2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], /** * Evolves the solution in time along the y axis (dimensional splitting) */ -__device__ void evolveG1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2], - float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1], - const int nx_, const int ny_, +template +__device__ void evolveG(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], + float G[block_height+1][block_width+1], const float dy_, const float dt_) { //Index of thread within block const int tx = threadIdx.x; const int ty = threadIdx.y; - //Index of cell within domain - const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1 - const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1; + const int i = tx + ghost_cells; //Skip local ghost cells, i.e., +1 + const int j = ty + ghost_cells; - if (ti > 0 && ti < nx_+1 && tj > 0 && tj < ny_+1) { - const int i = tx + 1; //Skip local ghost cells, i.e., +1 - const int j = ty + 1; - - Q[0][j][i] = Q[0][j][i] + (G[0][ty][tx] - G[0][ty+1][tx]) * dt_ / dy_; - Q[1][j][i] = Q[1][j][i] + (G[1][ty][tx] - G[1][ty+1][tx]) * dt_ / dy_; - Q[2][j][i] = Q[2][j][i] + (G[2][ty][tx] - G[2][ty+1][tx]) * dt_ / dy_; - } -} - - - - - - - -/** - * Evolves the solution in time along the y axis (dimensional splitting) - */ -__device__ void evolveG2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4], - float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1], - const int nx_, const int ny_, - const float dy_, const float dt_) { - //Index of thread within block - const int tx = threadIdx.x; - const int ty = threadIdx.y; - - //Index of cell within domain - const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2 - const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2; - - if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) { - const int i = tx + 2; //Skip local ghost cells, i.e., +2 - const int j = ty + 2; - - Q[0][j][i] = Q[0][j][i] + (G[0][ty][tx] - G[0][ty+1][tx]) * dt_ / dy_; - Q[1][j][i] = Q[1][j][i] + (G[1][ty][tx] - G[1][ty+1][tx]) * dt_ / dy_; - Q[2][j][i] = Q[2][j][i] + (G[2][ty][tx] - G[2][ty+1][tx]) * dt_ / dy_; - } + Q[j][i] = Q[j][i] + (G[ty][tx] - G[ty+1][tx]) * dt_ / dy_; }