diff --git a/GPUSimulators/Simulator.py b/GPUSimulators/Simulator.py
index d0c9f28..87319b4 100644
--- a/GPUSimulators/Simulator.py
+++ b/GPUSimulators/Simulator.py
@@ -83,9 +83,6 @@ class BaseSimulator:
         #Keep track of simulation time and number of timesteps
         self.t = 0.0
         self.nt = 0
-
-        #Log progress every n seconds during simulation
-        self.log_every = 5
         
                       
     def __str__(self):
@@ -102,31 +99,30 @@ class BaseSimulator:
     Requires that the stepEuler functionality is implemented in the subclasses
     """
     def simulateEuler(self, t_end):
-        with Common.Timer(self.__class__.__name__ + ".simulateEuler") as t:
-            # Compute number of timesteps to perform
-            n = int(t_end / self.dt + 1)
+        # Compute number of timesteps to perform
+        n = int(t_end / self.dt + 1)
 
-            next_print = self.log_every
+        printer = Common.ProgressPrinter(n)
+        
+        for i in range(0, n):
+            # Compute timestep for "this" iteration
+            local_dt = np.float32(min(self.dt, t_end-i*self.dt))
             
-            for i in range(0, n):
-                # Compute timestep for "this" iteration
-                local_dt = np.float32(min(self.dt, t_end-i*self.dt))
-                
-                # Stop if end reached (should not happen)
-                if (local_dt <= 0.0):
-                    break
-            
-                # Step with forward Euler 
-                self.stepEuler(local_dt)
+            # Stop if end reached (should not happen)
+            if (local_dt <= 0.0):
+                break
+        
+            # Step with forward Euler 
+            self.stepEuler(local_dt)
 
-                #Print info
-                if (t.elapsed() >= next_print):
-                    self.logger.info("%s simulated %d of %d steps (Euler)", self, i, n)
-                    next_print += self.log_every
-                    self.check()
+            #Print info
+            print_string = printer.getPrintString(i)
+            if (print_string):
+                self.logger.info("%s (Euler): %s", self, print_string)
+                self.check()
     
             
-        self.logger.info("%s simulated %f seconds to %f with %d steps (Euler)", self, t_end, self.t, n)
+        #self.logger.info("%s simulated %f seconds to %f with %d steps (Euler)", self, t_end, self.t, n)
         return self.t, n
         
     """
@@ -134,30 +130,28 @@ class BaseSimulator:
     Requires that the stepRK functionality is implemented in the subclasses
     """
     def simulateRK(self, t_end, order):
-        with Common.Timer(self.__class__.__name__ + ".simulateRK") as t:
-            # Compute number of timesteps to perform
-            n = int(t_end / self.dt + 1)
+        # Compute number of timesteps to perform
+        n = int(t_end / self.dt + 1)
+        
+        printer = Common.ProgressPrinter(n)
+        
+        for i in range(0, n):
+            # Compute timestep for "this" iteration
+            local_dt = np.float32(min(self.dt, t_end-i*self.dt))
+            
+            # Stop if end reached (should not happen)
+            if (local_dt <= 0.0):
+                break
+        
+            # Perform all the Runge-Kutta substeps
+            self.stepRK(local_dt, order)
 
-            next_print = self.log_every
-            
-            for i in range(0, n):
-                # Compute timestep for "this" iteration
-                local_dt = np.float32(min(self.dt, t_end-i*self.dt))
-                
-                # Stop if end reached (should not happen)
-                if (local_dt <= 0.0):
-                    break
-            
-                # Perform all the Runge-Kutta substeps
-                self.stepRK(local_dt, order)
-
-                #Print info
-                if (t.elapsed() >= next_print):
-                    self.logger.info("%s simulated %d of %d steps (RK2)", self, i, n)
-                    next_print += self.log_every
-                    self.check()
-            
-        self.logger.info("%s simulated %f seconds to %f with %d steps (RK2)", self, t_end, self.t, n)
+            #Print info
+            print_string = printer.getPrintString(i)
+            if (print_string):
+                self.logger.info("%s (RK2): %s", self, print_string)
+                self.check()
+    
         return self.t, n
         
     """
@@ -165,31 +159,29 @@ class BaseSimulator:
     Requires that the stepDimsplitX and stepDimsplitY functionality is implemented in the subclasses
     """
     def simulateDimsplit(self, t_end):
-        with Common.Timer(self.__class__.__name__ + ".simulateDimsplit") as t:
-            # Compute number of timesteps to perform
-            n = int(t_end / (2.0*self.dt) + 1)
-            
-            next_print = self.log_every
+        # Compute number of timesteps to perform
+        n = int(t_end / (2.0*self.dt) + 1)
+        
+        printer = Common.ProgressPrinter(n)
 
-            for i in range(0, n):
-                # Compute timestep for "this" iteration
-                local_dt = np.float32(0.5*min(2*self.dt, t_end-2*i*self.dt))
-                
-                # Stop if end reached (should not happen)
-                if (local_dt <= 0.0):
-                    break
-                
-                # Perform the dimensional split substeps
-                self.stepDimsplitXY(local_dt)
-                self.stepDimsplitYX(local_dt)
-
-                #Print info
-                if (t.elapsed() >= next_print):
-                    self.logger.info("%s simulated %d of %d steps (Dimsplit)", self, i, n)
-                    next_print += self.log_every
-                    self.check()
+        for i in range(0, n):
+            # Compute timestep for "this" iteration
+            local_dt = np.float32(0.5*min(2*self.dt, t_end-2*i*self.dt))
+            
+            # Stop if end reached (should not happen)
+            if (local_dt <= 0.0):
+                break
+            
+            # Perform the dimensional split substeps
+            self.stepDimsplitXY(local_dt)
+            self.stepDimsplitYX(local_dt)
+
+            #Print info
+            print_string = printer.getPrintString(i)
+            if (print_string):
+                self.logger.info("%s (Dimsplit): %s", self, print_string)
+                self.check()
             
-        self.logger.info("%s simulated %f seconds to %f with %d steps (Dimsplit)", self, t_end, self.t, 2*n)
         return self.t, 2*n
         
     
diff --git a/GPUSimulators/cuda/EulerCommon.h b/GPUSimulators/cuda/EulerCommon.h
index e34966a..02af22b 100644
--- a/GPUSimulators/cuda/EulerCommon.h
+++ b/GPUSimulators/cuda/EulerCommon.h
@@ -57,6 +57,59 @@ __device__ float4 F_func(const float4 Q, float P) {
 
 
 
+/**
+  * Harten-Lax-van Leer with contact discontinuity (Toro 2001, p 180)
+  */
+__device__ float4 HLL_flux(const float4 Q_l, const float4 Q_r, const float gamma) {
+    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;
+    
+    // Calculate pressures
+    const float P_l = pressure(Q_l, gamma);
+    const float P_r = pressure(Q_r, gamma);
+    
+    // Estimate the potential wave speeds
+    const float c_l = sqrt(gamma*P_l/Q_l.x);
+    const float c_r = sqrt(gamma*P_r/Q_r.x);
+    
+    // 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, P_l);
+    }
+    else if (S_r <= 0.0f) {
+        return F_func(Q_r, P_r);
+    }
+    //Or estimate flux in the star region
+    else {
+        const float4 F_l = F_func(Q_l, P_l);
+        const float4 F_r = F_func(Q_r, P_r);
+        const float4 flux = (S_r*F_l - S_l*F_r + S_r*S_l*(Q_r - Q_l)) / (S_r-S_l);
+        return flux;
+    }
+}
+
+
+
+
+
+
 
 /**
   * Central upwind flux function