Cleaned up include structure for kernels

SWESimulators/Common.py

@@ -8,6 +8,31 @@ import pycuda.compiler as cuda_compiler
 import pycuda.gpuarray
 import pycuda.driver as cuda
 
+
+
+"""
+Class which keeps track of time spent for a section of code
+"""
+class Timer(object):
+    def __init__(self, tag, verbose=True):
+        self.verbose = verbose
+        self.tag = tag
+        
+    def __enter__(self):
+        self.start = time.time()
+        return self
+    
+    def __exit__(self, *args):
+        self.end = time.time()
+        self.secs = self.end - self.start
+        self.msecs = self.secs * 1000 # millisecs
+        if self.verbose:
+            print("=> " + self.tag + ' %f ms' % self.msecs)
+            
+            
+            
+            
+
 """
 Class which keeps track of the CUDA context and some helper functions
 """
@@ -74,33 +99,76 @@ class CudaContext(object):
     """
     Reads a text file and creates an OpenCL kernel from that
     """
-    def get_kernel(self, kernel_filename, block_width, block_height):
+    def get_kernel(self, kernel_filename, block_width, block_height, include_dirs=[], verbose=False):
+        if (verbose):
+            print("Compiling " + kernel_filename)
+        
         # Generate a kernel ID for our cache
         module_path = os.path.dirname(os.path.realpath(__file__))
         
-        kernel_hash = ""
+        num_includes = 0
+        max_includes = 100
+        include_dirs = include_dirs + [module_path]
+        kernel_hash = kernel_filename + "_" + str(block_width) + "x" + str(block_height)
         
+        with Timer("compiler", verbose=False) as timer:
             # Loop over file and includes, and check if something has changed
-        files = [kernel_filename]
+            files = [os.path.join(module_path, kernel_filename)]
             while len(files):
-            filename = os.path.join(module_path, files.pop())
-            modified = os.path.getmtime(filename)
-            with open(filename, "r") as file:
-                file_str = file.read()
-                file_hash = filename + "_" + str(hash(file_str)) + ":" + str(modified) + "--"
-                includes = re.findall('^\W*#include\W+(.+?)\W*$', file_str, re.M)
-                files = files + includes #WARNING FIXME This will not work with circular includes
             
+                if (num_includes > max_includes):
+                    raise("Maximum number of includes reached - circular include in {:}?".format(kernel_filename))
+            
+                filename = files.pop()
+                
+                if (verbose):
+                    print("`-> Hashing " + filename)
+                    
+                modified = os.path.getmtime(filename)
+                    
+                # Open the file
+                with open(filename, "r") as file:
+                
+                    # Search for #inclue <something> and also hash the file
+                    file_str = file.read()
+                    file_hash = "\nfile=" + filename + ":hash=" + str(hash(file_str)) + ":modified=" + str(modified)
+                    includes = re.findall('^\W*#include\W+(.+?)\W*$', file_str, re.M)
+                    
+                # Loop over everything that looks like an include
+                for include_file in includes:
+                    
+                    #Search through include directories for the file
+                    file_path = os.path.dirname(filename)
+                    for include_path in [file_path] + include_dirs:
+                    
+                        # If we find it, add it to list of files to check
+                        temp_path = os.path.join(include_path, include_file)
+                        if (os.path.isfile(temp_path)):
+                            files = files + [temp_path]
+                            num_includes = num_includes + 1 #For circular includes...
+                            break
+                
+                #Create the "advanced" hash for each kernel
                 kernel_hash = kernel_hash + file_hash
+        if (verbose):
+            print("`-> Hashed in " + str(timer.secs) + " seconds")
         
         # Recompile kernel if file or includes have changed
         if (kernel_hash not in self.kernels.keys()):
+            if (verbose):
+                print("`-> Kernel not in hash => compiling " + kernel_filename)
+                
             #Create define string
             define_string = "#define block_width " + str(block_width) + "\n"
             define_string += "#define block_height " + str(block_height) + "\n\n"
             
             kernel_string = define_string + '#include "' + os.path.join(module_path, kernel_filename) + '"'
-            self.kernels[kernel_hash] = cuda_compiler.SourceModule(kernel_string, include_dirs=[module_path])
+            
+            with Timer("compiler", verbose=False) as timer:
+                self.kernels[kernel_hash] = cuda_compiler.SourceModule(kernel_string, include_dirs=include_dirs)
+            if (verbose):
+                print("`-> Compiled in " + str(timer.secs) + " seconds")
+            
             
         return self.kernels[kernel_hash]
     
@@ -113,21 +181,6 @@ class CudaContext(object):
         
         
         
-class Timer(object):
-    def __init__(self, tag, verbose=True):
-        self.verbose = verbose
-        self.tag = tag
-        
-    def __enter__(self):
-        self.start = time.time()
-        return self
-    
-    def __exit__(self, *args):
-        self.end = time.time()
-        self.secs = self.end - self.start
-        self.msecs = self.secs * 1000 # millisecs
-        if self.verbose:
-            print("=> " + self.tag + ' %f ms' % self.msecs)
         
         
 

SWESimulators/FORCE_kernel.cu

@@ -20,6 +20,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "fluxes/FirstOrderCentered.cu"
 
 
 /**

SWESimulators/HLL2_kernel.cu

@@ -19,6 +19,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "limiters.cu"
+#include "fluxes/HartenLaxVanLeer.cu"
 
 
 
@@ -64,7 +66,7 @@ void computeFluxF(float Q[3][block_height+4][block_width+4],
             const float3 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, g_) - F_func(Q_lr, g_));
 
             // Compute flux based on prediction
-            const float3 flux = HLLC_flux(Q_l_bar, Q_r_bar, g_);
+            const float3 flux = HLL_flux(Q_l_bar, Q_r_bar, g_);
             
             //Write to shared memory
             F[0][j][i] = flux.x;

SWESimulators/HLL_kernel.cu

@@ -1,5 +1,5 @@
 /*
-This OpenCL kernel implements the HLL flux
+This GPU kernel implements the HLL flux
 
 Copyright (C) 2016  SINTEF ICT
 
@@ -20,6 +20,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "fluxes/HartenLaxVanLeer.cu"
 
 
 

SWESimulators/KP07_dimsplit_kernel.cu

@@ -25,6 +25,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "limiters.cu"
+#include "fluxes/CentralUpwind.cu"
 
 
 __device__

SWESimulators/KP07_kernel.cu

@@ -25,6 +25,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "limiters.cu"
+#include "fluxes/CentralUpwind.cu"
 
 
 __device__

SWESimulators/LxF_kernel.cu

@@ -20,6 +20,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "fluxes/LaxFriedrichs.cu"
 
 
 /**

SWESimulators/WAF_kernel.cu

@@ -25,6 +25,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "common.cu"
+#include "fluxes/WeightedAverageFlux.cu"
 
 
 

SWESimulators/common.cu

@@ -464,79 +464,6 @@ __device__ void evolveG2(float Q[3][block_height+4][block_width+4],
 
 
 
-/**
-  * Reconstructs a slope using the minmod limiter based on three 
-  * consecutive values
-  */
-__device__ float minmodSlope(float left, float center, float right, float theta) {
-    const float backward = (center - left) * theta;
-    const float central = (right - left) * 0.5f;
-    const float forward = (right - center) * theta;
-    
-	return 0.25f
-		*copysign(1.0f, backward)
-		*(copysign(1.0f, backward) + copysign(1.0f, central))
-		*(copysign(1.0f, central) + copysign(1.0f, forward))
-		*min( min(fabs(backward), fabs(central)), fabs(forward) );
-}
-
-
-
-
-/**
-  * Reconstructs a minmod slope for a whole block along x
-  */
-__device__ void minmodSlopeX(float  Q[3][block_height+4][block_width+4],
-                  float Qx[3][block_height+2][block_width+2],
-                  const float theta_) {
-    //Index of thread within block
-    const int tx = get_local_id(0);
-    const int ty = get_local_id(1);
-    
-    //Reconstruct slopes along x axis
-    {
-        const int j = ty;
-        const int l = j + 2; //Skip ghost cells
-        for (int i=tx; i<block_width+2; i+=block_width) {
-            const int k = i + 1;
-            for (int p=0; p<3; ++p) {
-                Qx[p][j][i] = minmodSlope(Q[p][l][k-1], Q[p][l][k], Q[p][l][k+1], theta_);
-            }
-        }
-    }
-}
-
-
-/**
-  * Reconstructs a minmod slope for a whole block along y
-  */
-__device__ void minmodSlopeY(float  Q[3][block_height+4][block_width+4],
-                  float Qy[3][block_height+2][block_width+2],
-                  const float theta_) {
-    //Index of thread within block
-    const int tx = get_local_id(0);
-    const int ty = get_local_id(1);
-    
-    for (int j=ty; j<block_height+2; j+=block_height) {
-        const int l = j + 1;
-        {
-            const int i = tx;
-            const int k = i + 2; //Skip ghost cells
-            for (int p=0; p<3; ++p) {
-                Qy[p][j][i] = minmodSlope(Q[p][l-1][k], Q[p][l][k], Q[p][l+1][k], theta_);
-            }
-        }
-    }
-}
-
-
-
-
-
-
-
-
-
 __device__ float3 F_func(const float3 Q, const float g) {
     float3 F;
 
@@ -551,23 +478,9 @@ __device__ float3 F_func(const float3 Q, const float g) {
 
 
 
-/**
-  * Central upwind flux function
-  */
-__device__ float3 CentralUpwindFlux(const float3 Qm, float3 Qp, const float g) {
-    const float3 Fp = F_func(Qp, g);
-    const float up = Qp.y / Qp.x;   // hu / h
-    const float cp = sqrt(g*Qp.x); // sqrt(g*h)
 
-    const float3 Fm = F_func(Qm, g);
-    const float um = Qm.y / Qm.x;   // hu / h
-    const float cm = sqrt(g*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);
-}
 
 
 
@@ -578,397 +491,8 @@ __device__ float3 CentralUpwindFlux(const float3 Qm, float3 Qp, const float g) {
 
 
 
-/**
-  * Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 180)
-  */
-__device__ float3 HLL_flux(const float3 Q_l, const float3 Q_r, const float g_) {    
-    const float h_l = Q_l.x;
-    const float h_r = Q_r.x;
 
-    // Calculate velocities
-    const float u_l = Q_l.y / h_l;
-    const float u_r = Q_r.y / h_r;
 
-    // Estimate the potential wave speeds
-    const float c_l = sqrt(g_*h_l);
-    const float c_r = sqrt(g_*h_r);
-    
-    // Compute h in the "star region", h^dagger
-    const float h_dag = 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
-    
-    const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag / (h_l*h_l) ) );
-    const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag / (h_r*h_r) ) );
-    
-    const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
-    const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
-    
-    // Compute wave speed estimates
-    const float S_l = u_l - c_l*q_l;
-    const float S_r = u_r + c_r*q_r;
-    
-    //Upwind selection
-    if (S_l >= 0.0f) {
-        return F_func(Q_l, g_);
-    }
-    else if (S_r <= 0.0f) {
-        return F_func(Q_r, g_);
-    }
-    //Or estimate flux in the star region
-    else {
-        const float3 F_l = F_func(Q_l, g_);
-        const float3 F_r = F_func(Q_r, g_);
-        const float3 flux = (S_r*F_l - S_l*F_r + S_r*S_l*(Q_r - Q_l)) / (S_r-S_l);
-        return flux;
-    }
-}
-
-
-
-
-
-
-
-
-
-/**
-  * Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 181)
-  */
-__device__ float3 HLLC_flux(const float3 Q_l, const float3 Q_r, const float g_) {    
-    const float h_l = Q_l.x;
-    const float h_r = Q_r.x;
-    
-    // Calculate velocities
-    const float u_l = Q_l.y / h_l;
-    const float u_r = Q_r.y / h_r;
-    
-    // Estimate the potential wave speeds
-    const float c_l = sqrt(g_*h_l);
-    const float c_r = sqrt(g_*h_r);
-    
-    // Compute h in the "star region", h^dagger
-    const float h_dag = 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
-    
-    const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag / (h_l*h_l) ) );
-    const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag / (h_r*h_r) ) );
-    
-    const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
-    const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
-    
-    // Compute wave speed estimates
-    const float S_l = u_l - c_l*q_l;
-    const float S_r = u_r + c_r*q_r;
-    const float S_star = ( S_l*h_r*(u_r - S_r) - S_r*h_l*(u_l - S_l) ) / ( h_r*(u_r - S_r) - h_l*(u_l - S_l) );
-    
-    const float3 F_l = F_func(Q_l, g_);
-    const float3 F_r = F_func(Q_r, g_);
-    
-    //Upwind selection
-    if (S_l >= 0.0f) {
-        return F_l;
-    }
-    else if (S_r <= 0.0f) {
-        return F_r;
-    }
-    //Or estimate flux in the "left star" region
-    else if (S_l <= 0.0f && 0.0f <=S_star) {
-        const float v_l = Q_l.z / h_l;
-        const float3 Q_star_l = h_l * (S_l - u_l) / (S_l - S_star) * make_float3(1, S_star, v_l);
-        const float3 flux = F_l + S_l*(Q_star_l - Q_l);
-        return flux;
-    }
-    //Or estimate flux in the "righ star" region
-    else if (S_star <= 0.0f && 0.0f <=S_r) {
-        const float v_r = Q_r.z / h_r;
-        const float3 Q_star_r = h_r * (S_r - u_r) / (S_r - S_star) * make_float3(1, S_star, v_r);
-        const float3 flux = F_r + S_r*(Q_star_r - Q_r);
-        return flux;
-    }
-    else {
-        return make_float3(-99999.9f, -99999.9f, -99999.9f); //Something wrong here
-    }
-}
-
-
-
-/**
-  * Superbee flux limiter for WAF.
-  * Related to superbee limiter so that WAF_superbee(r, c) = 1 - (1-|c|)*superbee(r)
-  * @param r_ the ratio of upwind change (see Toro 2001, p. 203/204)
-  * @param c_ the courant number for wave k, dt*S_k/dx
-  */
-__device__ float WAF_superbee(float r_, float c_) {
-    // r <= 0.0
-    if (r_ <= 0.0f) { 
-        return 1.0f;
-    }
-    // 0.0 <= r <= 1/2
-    else if (r_ <= 0.5f) { 
-        return 1.0f - 2.0f*(1.0f - fabsf(c_))*r_;
-    }
-    // 1/2 <= r <= 1
-    else if (r_ <= 1.0f) {
-        return fabs(c_);
-    }
-    // 1 <= r <= 2
-    else  if (r_ <= 2.0f) {
-        return 1.0f - (1.0f - fabsf(c_))*r_;
-    }
-    // r >= 2
-    else {
-        return 2.0f*fabsf(c_) - 1.0f;
-    }
-}
-
-
-
-
-__device__ float WAF_albada(float r_, float c_) {
-    if (r_ <= 0.0f) {
-        return 1.0f;
-    }
-    else {
-        return 1.0f - (1.0f - fabsf(c_)) * r_ * (1.0f + r_) / (1.0f + r_*r_);
-    }
-}
-
-__device__ float WAF_minbee(float r_, float c_) {
-    r_ = fmaxf(-1.0f, fminf(2.0f, r_));
-    if (r_ <= 0.0f) {
-        return 1.0f;
-    }
-    if (r_ >= 0.0f && r_ <= 1.0f) {
-        return 1.0f - (1.0f - fabsf(c_)) * r_;
-    }
-    else {
-        return fabsf(c_);
-    }
-}
-
-__device__ float WAF_minmod(float r_, float c_) {
-    return 1.0f - (1.0f - fabsf(c_)) * fmaxf(0.0f, fminf(1.0f, r_));
-}
-
-
-__device__ float minmod(float r_) {
-    return fmaxf(0.0f, fminf(1.0f, r_));
-}
-
-__device__ float generalized_minmod(float r_, float theta_) {
-    return fmaxf(0.0f, fminf(theta_*r_, fminf( (1.0f + r_) / 2.0f, theta_)));
-}
-
-__device__ float superbee(float r_) {
-    return fmaxf(0.0f, fmaxf(fminf(2.0f*r_, 1.0f), fminf(r_, 2.0f)));
-}
-
-__device__ float vanAlbada1(float r_) {
-    return (r_*r_ + r_) / (r_*r_ + 1.0f);
-}
-
-__device__ float vanLeer(float r_) {
-    return (r_ + fabsf(r_)) / (1.0f + fabsf(r_));
-}
-
-__device__ float limiterToWAFLimiter(float r_, float c_) {
-    return 1.0f - (1.0f - fabsf(c_))*r_;
-}
-
-__device__ float desingularize(float x_, float eps_) {
-    return copysign(1.0f, x_)*fmaxf(fabsf(x_), fminf(x_*x_/(2.0f*eps_)+0.5f*eps_, eps_));
-}
-
-// Compute h in the "star region", h^dagger
-__device__ float computeHStar(float h_l, float h_r, float u_l, float u_r, float c_l, float c_r, float g_) {
-    //return 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
-    
-    const float h_tmp = 0.5f * (c_l + c_r) + 0.25f * (u_l - u_r);
-    return h_tmp*h_tmp / g_;
-}
-
-/**
-  * Weighted average flux (Toro 2001, p 200) for interface {i+1/2}
-  * @param r_ The flux limiter parameter (see Toro 2001, p. 203)
-  * @param Q_l2 Q_{i-1}
-  * @param Q_l1 Q_{i}
-  * @param Q_r1 Q_{i+1}
-  * @param Q_r2 Q_{i+2}
-  */
-__device__ float3 WAF_1D_flux(const float3 Q_l2, const float3 Q_l1, const float3 Q_r1, const float3 Q_r2, const float g_, const float dx_, const float dt_) {     
-    const float h_l = Q_l1.x;
-    const float h_r = Q_r1.x;
-    
-    const float h_l2 = Q_l2.x;
-    const float h_r2 = Q_r2.x;
-    
-    // Calculate velocities
-    const float u_l = Q_l1.y / h_l;
-    const float u_r = Q_r1.y / h_r;
-    
-    const float u_l2 = Q_l2.y / h_l2;
-    const float u_r2 = Q_r2.y / h_r2;
-    
-    const float v_l = Q_l1.z / h_l;
-    const float v_r = Q_r1.z / h_r;
-    
-    const float v_l2 = Q_l2.z / h_l2;
-    const float v_r2 = Q_r2.z / h_r2;
-    
-    // Estimate the potential wave speeds
-    const float c_l = sqrt(g_*h_l);
-    const float c_r = sqrt(g_*h_r);
-    
-    const float c_l2 = sqrt(g_*h_l2);
-    const float c_r2 = sqrt(g_*h_r2);
-    
-    // Compute h in the "star region", h^dagger
-    const float h_dag_l = computeHStar(h_l2,  h_l, u_l2,  u_l, c_l2,  c_l, g_);
-    const float h_dag   = computeHStar( h_l,  h_r,  u_l,  u_r,  c_l,  c_r, g_);
-    const float h_dag_r = computeHStar( h_r, h_r2,  u_r, u_r2,  c_r, c_r2, g_);
-    
-    const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag ) ) / h_l;
-    const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag ) ) / h_r;
-    
-    const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
-    const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
-    
-    // Compute wave speed estimates
-    const float S_l = u_l - c_l*q_l; 
-    const float S_r = u_r + c_r*q_r;
-    const float S_star = ( S_l*h_r*(u_r - S_r) - S_r*h_l*(u_l - S_l) ) / ( h_r*(u_r - S_r) - h_l*(u_l - S_l) );
-    
-    const float3 Q_star_l = h_l * (S_l - u_l) / (S_l - S_star) * make_float3(1.0, S_star, v_l);
-    const float3 Q_star_r = h_r * (S_r - u_r) / (S_r - S_star) * make_float3(1.0, S_star, v_r);
-    
-    // Estimate the fluxes in the four regions
-    const float3 F_1 = F_func(Q_l1, g_);
-    const float3 F_4 = F_func(Q_r1, g_);
-    
-    const float3 F_2 = F_1 + S_l*(Q_star_l - Q_l1);
-    const float3 F_3 = F_4 + S_r*(Q_star_r - Q_r1);
-    //const float3 F_2 = F_func(Q_star_l, g_);
-    //const float3 F_3 = F_func(Q_star_r, g_);
-    
-    // Compute the courant numbers for the waves
-    const float c_1 = S_l * dt_ / dx_;
-    const float c_2 = S_star * dt_ / dx_;
-    const float c_3 = S_r * dt_ / dx_;
-    
-    // Compute the "upwind change" vectors for the i-3/2 and i+3/2 interfaces
-    const float eps = 1.0e-6f;
-    const float r_1 = desingularize( (c_1 > 0.0f) ? (h_dag_l - h_l2) : (h_dag_r - h_r), eps) / desingularize((h_dag - h_l), eps);
-    const float r_2 = desingularize( (c_2 > 0.0f) ? (v_l - v_l2) : (v_r2 - v_r), eps ) / desingularize((v_r - v_l), eps);
-    const float r_3 = desingularize( (c_3 > 0.0f) ? (h_l - h_dag_l) : (h_r2 - h_dag_r), eps ) / desingularize((h_r - h_dag), eps);
-    
-    /*
-    const float eps1 = 1.0e-6f;
-    const float eps2 = eps1;//*1.0e-1f;
-    const float rh_m = desingularize((h_l - h_l2), eps1) / desingularize((h_r - h_l), eps2);
-    const float rh_p = desingularize((h_r2 - h_r), eps1) / desingularize((h_r - h_l), eps2);
-    const float rv_m = desingularize((v_l - v_l2), eps1) / desingularize((v_r - v_l), eps2);
-    const float rv_p = desingularize((v_r2 - v_r), eps1) / desingularize((v_r - v_l), eps2);
-    
-    // Compute the r parameters for the flux limiter
-    const float rh_1 = (c_1 > 0.0f) ? rh_m : rh_p; 
-    const float rv_1 = (c_1 > 0.0f) ? rv_m : rv_p; 
-    
-    const float rh_2 = (c_2 > 0.0f) ? rh_m : rh_p; 
-    const float rv_2 = (c_2 > 0.0f) ? rv_m : rv_p;
-    
-    const float rh_3 = (c_3 > 0.0f) ? rh_m : rh_p;
-    const float rv_3 = (c_3 > 0.0f) ? rv_m : rv_p;
-    
-    const float r_1 = rh_1;
-    const float r_2 = rv_2;
-    const float r_3 = rh_3;
-    */
-    
-    // Compute the limiter
-    // We use h for the nonlinear waves, and v for the middle shear wave 
-    const float A_1 = copysign(1.0f, c_1) * limiterToWAFLimiter(generalized_minmod(r_1, 1.9f), c_1);
-    const float A_2 = copysign(1.0f, c_2) * limiterToWAFLimiter(generalized_minmod(r_2, 1.9f), c_2); //Middle shear wave 
-    const float A_3 = copysign(1.0f, c_3) * limiterToWAFLimiter(generalized_minmod(r_3, 1.9f), c_3);
-    
-    //Average the fluxes
-    const float3 flux = 0.5f*( F_1 + F_4 )
-                      - 0.5f*( A_1 * (F_2 - F_1)
-                             + A_2 * (F_3 - F_2)
-                             + A_3 * (F_4 - F_3) );
-
-    return flux;
-}
-
-
-
-
-
-
-
-/**
-  * Lax-Friedrichs flux (Toro 2001, p 163)
-  */
-__device__ float3 LxF_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
-    const float3 F_l = F_func(Q_l, g_);
-    const float3 F_r = F_func(Q_r, g_);
-    
-    return 0.5f*(F_l + F_r) + (dx_/(2.0f*dt_))*(Q_l - Q_r);
-}
-
-
-
-/**
-  * Lax-Friedrichs extended to 2D
-  */
-__device__ float3 LxF_2D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
-    const float3 F_l = F_func(Q_l, g_);
-    const float3 F_r = F_func(Q_r, g_);
-    
-    //Note numerical diffusion for 2D here (0.25)
-    return 0.5f*(F_l + F_r) + (dx_/(4.0f*dt_))*(Q_l - Q_r);
-}
-
-
-
-
-/**
-  * Richtmeyer / Two-step Lax-Wendroff flux (Toro 2001, p 164)
-  */
-__device__ float3 LxW2_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
-    const float3 F_l = F_func(Q_l, g_);
-    const float3 F_r = F_func(Q_r, g_);
-    
-    const float3 Q_lw2 = 0.5f*(Q_l + Q_r) + (dt_/(2.0f*dx_))*(F_l - F_r);
-    
-    return F_func(Q_lw2, g_);
-}
-
-
-
-
-
-
-/**
-  * Godunovs centered scheme (Toro 2001, p 165)
-  */
-__device__ float3 GodC_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
-    const float3 F_l = F_func(Q_l, g_);
-    const float3 F_r = F_func(Q_r, g_);
-    
-    const float3 Q_godc = 0.5f*(Q_l + Q_r) + (dt_/dx_)*(F_l - F_r);
-    
-    return F_func(Q_godc, g_);
-}
-    
-
-    
-    
-/**
-  * First Ordered Centered (Toro 2001, p.163)
-  */
-__device__ float3 FORCE_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
-    const float3 F_lf = LxF_1D_flux(Q_l, Q_r, g_, dx_, dt_);
-    const float3 F_lw2 = LxW2_1D_flux(Q_l, Q_r, g_, dx_, dt_);
-    return 0.5f*(F_lf + F_lw2);
-}
 
 
 

SWESimulators/fluxes/CentralUpwind.cu

@@ -0,0 +1,39 @@
+/*
+This file implements the Central upwind flux
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+/**
+  * Central upwind flux function
+  */
+__device__ float3 CentralUpwindFlux(const float3 Qm, float3 Qp, const float g) {
+    const float3 Fp = F_func(Qp, g);
+    const float up = Qp.y / Qp.x;   // hu / h
+    const float cp = sqrt(g*Qp.x); // sqrt(g*h)
+
+    const float3 Fm = F_func(Qm, g);
+    const float um = Qm.y / Qm.x;   // hu / h
+    const float cm = sqrt(g*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);
+}
+

SWESimulators/fluxes/FirstOrderCentered.cu

@@ -0,0 +1,33 @@
+/*
+This file implements the First ORder CEntered flux
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "./LaxFriedrichs.cu"
+#include "./LaxWendroff.cu"
+
+
+    
+    
+/**
+  * First Ordered Centered (Toro 2001, p.163)
+  */
+__device__ float3 FORCE_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
+    const float3 F_lf = LxF_1D_flux(Q_l, Q_r, g_, dx_, dt_);
+    const float3 F_lw2 = LxW2_1D_flux(Q_l, Q_r, g_, dx_, dt_);
+    return 0.5f*(F_lf + F_lw2);
+}

SWESimulators/fluxes/Godunov.cu

@@ -0,0 +1,36 @@
+/*
+This file implements the Godunov flux
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+
+/**
+  * Godunovs centered scheme (Toro 2001, p 165)
+  */
+__device__ float3 GodC_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
+    const float3 F_l = F_func(Q_l, g_);
+    const float3 F_r = F_func(Q_r, g_);
+    
+    const float3 Q_godc = 0.5f*(Q_l + Q_r) + (dt_/dx_)*(F_l - F_r);
+    
+    return F_func(Q_godc, g_);
+}
+    
+
+

SWESimulators/fluxes/HartenLaxVanLeer.cu

@@ -0,0 +1,63 @@
+/*
+This file implements the Harten-Lax-van Leer flux
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+
+/**
+  * Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 180)
+  */
+__device__ float3 HLL_flux(const float3 Q_l, const float3 Q_r, const float g_) {    
+    const float h_l = Q_l.x;
+    const float h_r = Q_r.x;
+    
+    // Calculate velocities
+    const float u_l = Q_l.y / h_l;
+    const float u_r = Q_r.y / h_r;
+    
+    // Estimate the potential wave speeds
+    const float c_l = sqrt(g_*h_l);
+    const float c_r = sqrt(g_*h_r);
+    
+    // Compute h in the "star region", h^dagger
+    const float h_dag = 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
+    
+    const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag / (h_l*h_l) ) );
+    const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag / (h_r*h_r) ) );
+    
+    const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
+    const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
+    
+    // Compute wave speed estimates
+    const float S_l = u_l - c_l*q_l;
+    const float S_r = u_r + c_r*q_r;
+    
+    //Upwind selection
+    if (S_l >= 0.0f) {
+        return F_func(Q_l, g_);
+    }
+    else if (S_r <= 0.0f) {
+        return F_func(Q_r, g_);
+    }
+    //Or estimate flux in the star region
+    else {
+        const float3 F_l = F_func(Q_l, g_);
+        const float3 F_r = F_func(Q_r, g_);
+        const float3 flux = (S_r*F_l - S_l*F_r + S_r*S_l*(Q_r - Q_l)) / (S_r-S_l);
+        return flux;
+    }
+}

SWESimulators/fluxes/HartenLaxVanLeerContact.cu

@@ -0,0 +1,80 @@
+/*
+This file implements the Harten-Lax-van Leer flux with contact discontinuity
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+
+/**
+  * Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 181)
+  */
+__device__ float3 HLLC_flux(const float3 Q_l, const float3 Q_r, const float g_) {    
+    const float h_l = Q_l.x;
+    const float h_r = Q_r.x;
+    
+    // Calculate velocities
+    const float u_l = Q_l.y / h_l;
+    const float u_r = Q_r.y / h_r;
+    
+    // Estimate the potential wave speeds
+    const float c_l = sqrt(g_*h_l);
+    const float c_r = sqrt(g_*h_r);
+    
+    // Compute h in the "star region", h^dagger
+    const float h_dag = 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
+    
+    const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag / (h_l*h_l) ) );
+    const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag / (h_r*h_r) ) );
+    
+    const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
+    const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
+    
+    // Compute wave speed estimates
+    const float S_l = u_l - c_l*q_l;
+    const float S_r = u_r + c_r*q_r;
+    const float S_star = ( S_l*h_r*(u_r - S_r) - S_r*h_l*(u_l - S_l) ) / ( h_r*(u_r - S_r) - h_l*(u_l - S_l) );
+    
+    const float3 F_l = F_func(Q_l, g_);
+    const float3 F_r = F_func(Q_r, g_);
+    
+    //Upwind selection
+    if (S_l >= 0.0f) {
+        return F_l;
+    }
+    else if (S_r <= 0.0f) {
+        return F_r;
+    }
+    //Or estimate flux in the "left star" region
+    else if (S_l <= 0.0f && 0.0f <=S_star) {
+        const float v_l = Q_l.z / h_l;
+        const float3 Q_star_l = h_l * (S_l - u_l) / (S_l - S_star) * make_float3(1, S_star, v_l);
+        const float3 flux = F_l + S_l*(Q_star_l - Q_l);
+        return flux;
+    }
+    //Or estimate flux in the "righ star" region
+    else if (S_star <= 0.0f && 0.0f <=S_r) {
+        const float v_r = Q_r.z / h_r;
+        const float3 Q_star_r = h_r * (S_r - u_r) / (S_r - S_star) * make_float3(1, S_star, v_r);
+        const float3 flux = F_r + S_r*(Q_star_r - Q_r);
+        return flux;
+    }
+    else {
+        return make_float3(-99999.9f, -99999.9f, -99999.9f); //Something wrong here
+    }
+}
+

SWESimulators/fluxes/LaxFriedrichs.cu

@@ -0,0 +1,24 @@
+
+
+/**
+  * Lax-Friedrichs flux (Toro 2001, p 163)
+  */
+__device__ float3 LxF_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
+    const float3 F_l = F_func(Q_l, g_);
+    const float3 F_r = F_func(Q_r, g_);
+    
+    return 0.5f*(F_l + F_r) + (dx_/(2.0f*dt_))*(Q_l - Q_r);
+}
+
+
+
+/**
+  * Lax-Friedrichs extended to 2D
+  */
+__device__ float3 LxF_2D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
+    const float3 F_l = F_func(Q_l, g_);
+    const float3 F_r = F_func(Q_r, g_);
+    
+    //Note numerical diffusion for 2D here (0.25)
+    return 0.5f*(F_l + F_r) + (dx_/(4.0f*dt_))*(Q_l - Q_r);
+}

SWESimulators/fluxes/LaxWendroff.cu

@@ -0,0 +1,33 @@
+/*
+This file implements the Lax-Wendroff flux
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+/**
+  * Richtmeyer / Two-step Lax-Wendroff flux (Toro 2001, p 164)
+  */
+__device__ float3 LxW2_1D_flux(const float3 Q_l, const float3 Q_r, const float g_, const float dx_, const float dt_) {
+    const float3 F_l = F_func(Q_l, g_);
+    const float3 F_r = F_func(Q_r, g_);
+    
+    const float3 Q_lw2 = 0.5f*(Q_l + Q_r) + (dt_/(2.0f*dx_))*(F_l - F_r);
+    
+    return F_func(Q_lw2, g_);
+}
+

SWESimulators/fluxes/WeightedAverageFlux.cu

@@ -0,0 +1,187 @@
+/*
+This file implements the Weighted Average Flux
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "limiters.cu"
+
+/**
+  * Superbee flux limiter for WAF.
+  * Related to superbee limiter so that WAF_superbee(r, c) = 1 - (1-|c|)*superbee(r)
+  * @param r_ the ratio of upwind change (see Toro 2001, p. 203/204)
+  * @param c_ the courant number for wave k, dt*S_k/dx
+  */
+__device__ float WAF_superbee(float r_, float c_) {
+    // r <= 0.0
+    if (r_ <= 0.0f) { 
+        return 1.0f;
+    }
+    // 0.0 <= r <= 1/2
+    else if (r_ <= 0.5f) { 
+        return 1.0f - 2.0f*(1.0f - fabsf(c_))*r_;
+    }
+    // 1/2 <= r <= 1
+    else if (r_ <= 1.0f) {
+        return fabs(c_);
+    }
+    // 1 <= r <= 2
+    else  if (r_ <= 2.0f) {
+        return 1.0f - (1.0f - fabsf(c_))*r_;
+    }
+    // r >= 2
+    else {
+        return 2.0f*fabsf(c_) - 1.0f;
+    }
+}
+
+
+
+
+__device__ float WAF_albada(float r_, float c_) {
+    if (r_ <= 0.0f) {
+        return 1.0f;
+    }
+    else {
+        return 1.0f - (1.0f - fabsf(c_)) * r_ * (1.0f + r_) / (1.0f + r_*r_);
+    }
+}
+
+__device__ float WAF_minbee(float r_, float c_) {
+    r_ = fmaxf(-1.0f, fminf(2.0f, r_));
+    if (r_ <= 0.0f) {
+        return 1.0f;
+    }
+    if (r_ >= 0.0f && r_ <= 1.0f) {
+        return 1.0f - (1.0f - fabsf(c_)) * r_;
+    }
+    else {
+        return fabsf(c_);
+    }
+}
+
+__device__ float WAF_minmod(float r_, float c_) {
+    return 1.0f - (1.0f - fabsf(c_)) * fmaxf(0.0f, fminf(1.0f, r_));
+}
+
+
+
+__device__ float limiterToWAFLimiter(float r_, float c_) {
+    return 1.0f - (1.0f - fabsf(c_))*r_;
+}
+
+__device__ float desingularize(float x_, float eps_) {
+    return copysign(1.0f, x_)*fmaxf(fabsf(x_), fminf(x_*x_/(2.0f*eps_)+0.5f*eps_, eps_));
+}
+
+// Compute h in the "star region", h^dagger
+__device__ __inline__ float computeHStar(float h_l, float h_r, float u_l, float u_r, float c_l, float c_r, float g_) {
+    
+    //This estimate for the h* gives rise to spurious oscillations. 
+    //return 0.5f * (h_l+h_r) - 0.25f * (u_r-u_l)*(h_l+h_r)/(c_l+c_r);
+    
+    const float h_tmp = 0.5f * (c_l + c_r) + 0.25f * (u_l - u_r);
+    return h_tmp*h_tmp / g_;
+}
+
+/**
+  * Weighted average flux (Toro 2001, p 200) for interface {i+1/2}
+  * @param r_ The flux limiter parameter (see Toro 2001, p. 203)
+  * @param Q_l2 Q_{i-1}
+  * @param Q_l1 Q_{i}
+  * @param Q_r1 Q_{i+1}
+  * @param Q_r2 Q_{i+2}
+  */
+__device__ float3 WAF_1D_flux(const float3 Q_l2, const float3 Q_l1, const float3 Q_r1, const float3 Q_r2, const float g_, const float dx_, const float dt_) {     
+    const float h_l = Q_l1.x;
+    const float h_r = Q_r1.x;
+    
+    const float h_l2 = Q_l2.x;
+    const float h_r2 = Q_r2.x;
+    
+    // Calculate velocities
+    const float u_l = Q_l1.y / h_l;
+    const float u_r = Q_r1.y / h_r;
+    
+    const float u_l2 = Q_l2.y / h_l2;
+    const float u_r2 = Q_r2.y / h_r2;
+    
+    const float v_l = Q_l1.z / h_l;
+    const float v_r = Q_r1.z / h_r;
+    
+    const float v_l2 = Q_l2.z / h_l2;
+    const float v_r2 = Q_r2.z / h_r2;
+    
+    // Estimate the potential wave speeds
+    const float c_l = sqrt(g_*h_l);
+    const float c_r = sqrt(g_*h_r);
+    
+    const float c_l2 = sqrt(g_*h_l2);
+    const float c_r2 = sqrt(g_*h_r2);
+    
+    // Compute h in the "star region", h^dagger
+    const float h_dag_l = computeHStar(h_l2,  h_l, u_l2,  u_l, c_l2,  c_l, g_);
+    const float h_dag   = computeHStar( h_l,  h_r,  u_l,  u_r,  c_l,  c_r, g_);
+    const float h_dag_r = computeHStar( h_r, h_r2,  u_r, u_r2,  c_r, c_r2, g_);
+    
+    const float q_l_tmp = sqrt(0.5f * ( (h_dag+h_l)*h_dag ) ) / h_l;
+    const float q_r_tmp = sqrt(0.5f * ( (h_dag+h_r)*h_dag ) ) / h_r;
+    
+    const float q_l = (h_dag > h_l) ? q_l_tmp : 1.0f;
+    const float q_r = (h_dag > h_r) ? q_r_tmp : 1.0f;
+    
+    // Compute wave speed estimates
+    const float S_l = u_l - c_l*q_l; 
+    const float S_r = u_r + c_r*q_r;
+    const float S_star = ( S_l*h_r*(u_r - S_r) - S_r*h_l*(u_l - S_l) ) / ( h_r*(u_r - S_r) - h_l*(u_l - S_l) );
+    
+    const float3 Q_star_l = h_l * (S_l - u_l) / (S_l - S_star) * make_float3(1.0, S_star, v_l);
+    const float3 Q_star_r = h_r * (S_r - u_r) / (S_r - S_star) * make_float3(1.0, S_star, v_r);
+    
+    // Estimate the fluxes in the four regions
+    const float3 F_1 = F_func(Q_l1, g_);
+    const float3 F_4 = F_func(Q_r1, g_);
+    
+    const float3 F_2 = F_1 + S_l*(Q_star_l - Q_l1);
+    const float3 F_3 = F_4 + S_r*(Q_star_r - Q_r1);
+    //const float3 F_2 = F_func(Q_star_l, g_);
+    //const float3 F_3 = F_func(Q_star_r, g_);
+    
+    // Compute the courant numbers for the waves
+    const float c_1 = S_l * dt_ / dx_;
+    const float c_2 = S_star * dt_ / dx_;
+    const float c_3 = S_r * dt_ / dx_;
+    
+    // Compute the "upwind change" vectors for the i-3/2 and i+3/2 interfaces
+    const float eps = 1.0e-6f;
+    const float r_1 = desingularize( (c_1 > 0.0f) ? (h_dag_l - h_l2) : (h_dag_r - h_r), eps) / desingularize((h_dag - h_l), eps);
+    const float r_2 = desingularize( (c_2 > 0.0f) ? (v_l - v_l2) : (v_r2 - v_r), eps ) / desingularize((v_r - v_l), eps);
+    const float r_3 = desingularize( (c_3 > 0.0f) ? (h_l - h_dag_l) : (h_r2 - h_dag_r), eps ) / desingularize((h_r - h_dag), eps);
+        
+    // Compute the limiter
+    // We use h for the nonlinear waves, and v for the middle shear wave 
+    const float A_1 = copysign(1.0f, c_1) * limiterToWAFLimiter(generalized_minmod(r_1, 1.9f), c_1);
+    const float A_2 = copysign(1.0f, c_2) * limiterToWAFLimiter(generalized_minmod(r_2, 1.9f), c_2); 
+    const float A_3 = copysign(1.0f, c_3) * limiterToWAFLimiter(generalized_minmod(r_3, 1.9f), c_3);
+    
+    //Average the fluxes
+    const float3 flux = 0.5f*( F_1 + F_4 )
+                      - 0.5f*( A_1 * (F_2 - F_1)
+                             + A_2 * (F_3 - F_2)
+                             + A_3 * (F_4 - F_3) );
+
+    return flux;
+}

SWESimulators/limiters.cu

@@ -0,0 +1,129 @@
+/*
+This file implements different flux and slope limiters
+
+Copyright (C) 2016, 2017, 2018 SINTEF ICT
+
+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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+
+
+
+
+
+/**
+  * Reconstructs a slope using the generalized minmod limiter based on three 
+  * consecutive values
+  */
+__device__ __inline__ float minmodSlope(float left, float center, float right, float theta) {
+    const float backward = (center - left) * theta;
+    const float central = (right - left) * 0.5f;
+    const float forward = (right - center) * theta;
+    
+	return 0.25f
+		*copysign(1.0f, backward)
+		*(copysign(1.0f, backward) + copysign(1.0f, central))
+		*(copysign(1.0f, central) + copysign(1.0f, forward))
+		*min( min(fabs(backward), fabs(central)), fabs(forward) );
+}
+
+
+
+
+/**
+  * Reconstructs a minmod slope for a whole block along the abscissa
+  */
+__device__ void minmodSlopeX(float  Q[3][block_height+4][block_width+4],
+                  float Qx[3][block_height+2][block_width+2],
+                  const float theta_) {
+    //Index of thread within block
+    const int tx = get_local_id(0);
+    const int ty = get_local_id(1);
+    
+    //Reconstruct slopes along x axis
+    {
+        const int j = ty;
+        const int l = j + 2; //Skip ghost cells
+        for (int i=tx; i<block_width+2; i+=block_width) {
+            const int k = i + 1;
+            for (int p=0; p<3; ++p) {
+                Qx[p][j][i] = minmodSlope(Q[p][l][k-1], Q[p][l][k], Q[p][l][k+1], theta_);
+            }
+        }
+    }
+}
+
+
+/**
+  * Reconstructs a minmod slope for a whole block along the ordinate
+  */
+__device__ void minmodSlopeY(float  Q[3][block_height+4][block_width+4],
+                  float Qy[3][block_height+2][block_width+2],
+                  const float theta_) {
+    //Index of thread within block
+    const int tx = get_local_id(0);
+    const int ty = get_local_id(1);
+    
+    for (int j=ty; j<block_height+2; j+=block_height) {
+        const int l = j + 1;
+        {
+            const int i = tx;
+            const int k = i + 2; //Skip ghost cells
+            for (int p=0; p<3; ++p) {
+                Qy[p][j][i] = minmodSlope(Q[p][l-1][k], Q[p][l][k], Q[p][l+1][k], theta_);
+            }
+        }
+    }
+}
+
+
+
+
+__device__ float monotonized_central(float r_) {
+    return fmaxf(0.0f, fminf(2.0f, fminf(2.0f*r_, 0.5f*(1.0f+r_))));
+}
+
+__device__ float osher(float r_, float beta_) {
+    return fmaxf(0.0f, fminf(beta_, r_));
+}
+
+__device__ float sweby(float r_, float beta_) {
+    return fmaxf(0.0f, fmaxf(fminf(r_, beta_), fminf(beta_*r_, 1.0f)));
+}
+
+__device__ float minmod(float r_) {
+    return fmaxf(0.0f, fminf(1.0f, r_));
+}
+
+__device__ float generalized_minmod(float r_, float theta_) {
+    return fmaxf(0.0f, fminf(theta_*r_, fminf( (1.0f + r_) / 2.0f, theta_)));
+}
+
+__device__ float superbee(float r_) {
+    return fmaxf(0.0f, fmaxf(fminf(2.0f*r_, 1.0f), fminf(r_, 2.0f)));
+}
+
+__device__ float vanAlbada1(float r_) {
+    return (r_*r_ + r_) / (r_*r_ + 1.0f);
+}
+
+__device__ float vanAlbada2(float r_) {
+    return 2.0f*r_ / (r_*r_* + 1.0f);
+}
+
+__device__ float vanLeer(float r_) {
+    return (r_ + fabsf(r_)) / (1.0f + fabsf(r_));
+}