diff --git a/GPUSimulators/helpers/InitialConditions.py b/GPUSimulators/helpers/InitialConditions.py
new file mode 100644
index 0000000..724823c
--- /dev/null
+++ b/GPUSimulators/helpers/InitialConditions.py
@@ -0,0 +1,261 @@
+# -*- coding: utf-8 -*-
+
+"""
+This python module implements Cuda context handling
+
+Copyright (C) 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/>.
+"""
+
+
+from GPUSimulators.Simulator import BoundaryCondition
+import numpy as np
+
+def genShockBubble(nx, ny, gamma):
+    """
+    Generate Shock-bubble interaction case for the Euler equations
+    """
+    
+    width = 4.0
+    height = 1.0
+    dx = width / float(nx)
+    dy = height / float(ny)
+    g = 0.0
+
+
+    rho = np.ones((ny, nx), dtype=np.float32)
+    u = np.zeros((ny, nx), dtype=np.float32)
+    v = np.zeros((ny, nx), dtype=np.float32)
+    E = np.zeros((ny, nx), dtype=np.float32)
+    p = np.ones((ny, nx), dtype=np.float32)
+
+    x_center = 0.5
+    y_center = 0.5
+    x = np.linspace(0.5*dx, nx*dx-0.5*dx, nx, dtype=np.float32) - x_center
+    y = np.linspace(0.5*dy, ny*dy-0.5*dy, ny, dtype=np.float32) - y_center
+    xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
+       
+    #Bubble
+    radius = 0.25
+    bubble = np.sqrt(xv**2+yv**2) <= radius
+    rho = np.where(bubble, 0.1, rho)
+    
+    #Left boundary
+    left = (xv - xv.min() < 0.1)
+    rho = np.where(left, 3.81250, rho)
+    u = np.where(left, 2.57669, u)
+    
+    #Energy
+    p = np.where(left, 10.0, p)
+    E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)
+    
+    #Estimate dt
+    scale = 0.45
+    max_rho_estimate = 5.0
+    max_u_estimate = 5.0
+    dx = width/nx
+    dy = height/ny
+    dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))
+
+    bc = BoundaryCondition({
+        'north': BoundaryCondition.Type.Reflective,
+        'south': BoundaryCondition.Type.Reflective,
+        'east': BoundaryCondition.Type.Periodic,
+        'west': BoundaryCondition.Type.Periodic
+    })
+    
+    #Construct simulator
+    arguments = {
+        'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
+        'nx': nx, 'ny': ny,
+        'dx': dx, 'dy': dy, 'dt': dt,
+        'g': g,
+        'gamma': gamma,
+        'boundary_conditions': bc
+    } 
+    return arguments
+
+    
+    
+    
+    
+    
+    
+def genKelvinHelmholtz(nx, ny, gamma, roughness=0.125):
+    """
+    Roughness parameter in (0, 1.0] determines how "squiggly" 
+    the interface betweeen the zones is
+    """
+    
+    def genZones(nx, ny, n):
+        """
+        Generates the zones of the two fluids of K-H
+        """
+        zone = np.zeros((ny, nx), dtype=np.int32)
+
+        dx = 1.0 / nx
+        dy = 1.0 / ny
+
+        def genSmoothRandom(nx, n):
+            assert (n <= nx), "Number of generated points nx must be larger than n"
+            
+            if n == nx:
+                return np.random.random(nx)-0.5
+            else:
+                from scipy.interpolate import interp1d
+
+                #Control points and interpolator
+                xp = np.linspace(0.0, 1.0, n)
+                yp = np.random.random(n) - 0.5
+                f = interp1d(xp, yp, kind='cubic')
+
+                #Interpolation points
+                x = np.linspace(0.0, 1.0, nx)
+                return f(x)
+
+
+
+        x = np.linspace(0, 1, nx)
+        y = np.linspace(0, 1, ny)
+
+        _, y = np.meshgrid(x, y)
+
+        #print(y+a[0])
+
+        a = genSmoothRandom(nx, n)*dy
+        zone = np.where(y > 0.25+a, zone, 1)
+
+        a = genSmoothRandom(nx, n)*dy
+        zone = np.where(y < 0.75+a, zone, 1)
+        
+        return zone
+        
+    width = 2.0
+    height = 1.0
+    dx = width / float(nx)
+    dy = height / float(ny)
+    g = 0.0
+    gamma = 1.4
+
+    rho = np.empty((ny, nx), dtype=np.float32)
+    u = np.empty((ny, nx), dtype=np.float32)
+    v = np.zeros((ny, nx), dtype=np.float32)
+    p = 2.5*np.ones((ny, nx), dtype=np.float32)
+
+    #Generate the different zones    
+    zones = genZones(nx, ny, max(1, min(nx, int(nx*roughness))))
+    
+    #Zone 0
+    zone0 = zones == 0
+    rho = np.where(zone0, 1.0, rho)
+    u = np.where(zone0, 0.5, u)
+    
+    #Zone 1
+    zone1 = zones == 1
+    rho = np.where(zone1, 2.0, rho)
+    u = np.where(zone1, -0.5, u)
+    
+    E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)
+    
+    #Estimate dt
+    scale = 0.9
+    max_rho_estimate = 3.0
+    max_u_estimate = 2.0
+    dx = width/nx
+    dy = height/ny
+    dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))
+    
+    
+    bc = BoundaryCondition({
+        'north': BoundaryCondition.Type.Periodic,
+        'south': BoundaryCondition.Type.Periodic,
+        'east': BoundaryCondition.Type.Periodic,
+        'west': BoundaryCondition.Type.Periodic
+    })
+    
+    #Construct simulator
+    arguments = {
+        'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
+        'nx': nx, 'ny': ny,
+        'dx': dx, 'dy': dy, 'dt': dt,
+        'g': g,
+        'gamma': gamma,
+        'boundary_conditions': bc
+    } 
+    
+    return arguments
+    
+    
+    
+def genRayleighTaylor(nx, ny, gamma, version=0):
+    """
+    Generates Rayleigh-Taylor instability case
+    """
+    width = 0.5
+    height = 1.5
+    dx = width / float(nx)
+    dy = height / float(ny)
+    g = 0.1
+
+    rho = np.zeros((ny, nx), dtype=np.float32)
+    u = np.zeros((ny, nx), dtype=np.float32)
+    v = np.zeros((ny, nx), dtype=np.float32)
+    p = np.zeros((ny, nx), dtype=np.float32)
+    
+    x = np.linspace(0.5*dx, nx*dx-0.5*dx, nx, dtype=np.float32)-width*0.5
+    y = np.linspace(0.5*dy, ny*dy-0.5*dy, ny, dtype=np.float32)-height*0.5
+    xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
+    
+    #This gives a squigly interfact
+    if (version == 0):
+        y_threshold = 0.01*np.cos(2*np.pi*np.abs(x)/0.5)
+        rho = np.where(yv <= y_threshold, 1.0, rho)
+        rho = np.where(yv > y_threshold, 2.0, rho)
+    elif (version == 1):
+        rho = np.where(yv <= 0.0, 1.0, rho)
+        rho = np.where(yv > 0.0, 2.0, rho)
+        v = 0.01*(1.0 + np.cos(2*np.pi*xv/0.5))/4
+    else:
+        assert False, "Invalid version"
+    
+    p = 2.5 - rho*g*yv
+    E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)
+    
+    #Estimate dt
+    scale = 0.9
+    max_rho_estimate = 3.0
+    max_u_estimate = 1.0
+    dx = width/nx
+    dy = height/ny
+    dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))
+    
+    bc = BoundaryCondition({
+        'north': BoundaryCondition.Type.Reflective,
+        'south': BoundaryCondition.Type.Reflective,
+        'east': BoundaryCondition.Type.Reflective,
+        'west': BoundaryCondition.Type.Reflective
+    })
+    
+    #Construct simulator
+    arguments = {
+        'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
+        'nx': nx, 'ny': ny,
+        'dx': dx, 'dy': dy, 'dt': dt,
+        'g': g,
+        'gamma': gamma,
+        'boundary_conditions': bc
+    } 
+
+    return arguments
\ No newline at end of file
diff --git a/GPUSimulators/helpers/Visualization.py b/GPUSimulators/helpers/Visualization.py
new file mode 100644
index 0000000..8314952
--- /dev/null
+++ b/GPUSimulators/helpers/Visualization.py
@@ -0,0 +1,61 @@
+# -*- coding: utf-8 -*-
+
+"""
+This python module implements Cuda context handling
+
+Copyright (C) 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/>.
+"""
+
+
+
+import numpy as np
+
+from matplotlib.colors import Normalize
+
+
+
+def genSchlieren(rho):
+    #Compute length of z-component of normalized gradient vector 
+    normal = np.gradient(rho) #[x, y, 1]
+    length = 1.0 / np.sqrt(normal[0]**2 + normal[1]**2 + 1.0)
+    schlieren = np.power(length, 128)
+    return schlieren
+
+
+def genVorticity(rho, rho_u, rho_v):
+    u = rho_u / rho
+    v = rho_v / rho
+    u = np.sqrt(u**2 + v**2)
+    u_max = u.max()
+    
+    du_dy, _ = np.gradient(u)
+    _, dv_dx = np.gradient(v)
+    
+    #Length of curl
+    curl = dv_dx - du_dy
+    return curl
+
+
+def genColors(rho, rho_u, rho_v, cmap, vmax, vmin):
+    schlieren = genSchlieren(rho)
+    curl = genVorticity(rho, rho_u, rho_v)
+
+    colors = Normalize(vmin, vmax, clip=True)(curl)
+    colors = cmap(colors)
+    for k in range(3):
+        colors[:,:,k] = colors[:,:,k]*schlieren
+
+    return colors
\ No newline at end of file