Bugfix KP07 and refactoring

GPUSimulators/EE2D_KP07_dimsplit.py

@@ -102,8 +102,8 @@ class EE2D_KP07_dimsplit (BaseSimulator):
                         2, 2, 
                         [None, None, None, None])
         self.cfl_data = gpuarray.GPUArray(self.grid_size, dtype=np.float32)
-        dt_x = np.min(self.dx / (np.abs(hu0/h0) + np.sqrt(gamma*h0)))
-        dt_y = np.min(self.dy / (np.abs(hv0/h0) + np.sqrt(gamma*h0)))
+        dt_x = np.min(self.dx / (np.abs(rho_u/rho) + np.sqrt(gamma*rho)))
+        dt_y = np.min(self.dy / (np.abs(rho_v/rho) + np.sqrt(gamma*rho)))
         dt = min(dt_x, dt_y)
         self.cfl_data.fill(dt, stream=self.stream)
                         

GPUSimulators/LxF.py

@@ -117,6 +117,10 @@ class LxF (Simulator.BaseSimulator):
     def download(self):
         return self.u0.download(self.stream)
 
+    def check(self):
+        self.u0.check()
+        self.u1.check()
+        
     def computeDt(self):
         max_dt = gpuarray.min(self.cfl_data, stream=self.stream).get();
         return max_dt*0.5

GPUSimulators/Simulator.py

@@ -159,7 +159,7 @@ class BaseSimulator(object):
         return "{:s} [{:d}x{:d}]".format(self.__class__.__name__, self.nx, self.ny)
 
 
-    def simulate(self, t):
+    def simulate(self, t, dt=None):
         """ 
         Function which simulates t_end seconds using the step function
         Requires that the step() function is implemented in the subclasses
@@ -167,27 +167,31 @@ class BaseSimulator(object):
 
         printer = Common.ProgressPrinter(t)
         
-        t_end = self.simTime() + t
+        t_start = self.simTime()
+        t_end = t_start + t
         
-        dt = None
+        local_dt = dt
+        
+        if (local_dt == None):
+            local_dt = self.computeDt()
         
         while(self.simTime() < t_end):
-            if (self.simSteps() % 100 == 0):
-                dt = self.computeDt()*self.cfl_scale
+            if (dt == None and self.simSteps() % 100 == 0):
+                local_dt = self.computeDt()
         
             # Compute timestep for "this" iteration (i.e., shorten last timestep)
-            dt = np.float32(min(dt, t_end-self.simTime()))
+            local_dt = np.float32(min(local_dt*self.cfl_scale, t_end-self.simTime()))
 
             # Stop if end reached (should not happen)
-            if (dt <= 0.0):
+            if (local_dt <= 0.0):
                 self.logger.warning("Timestep size {:d} is less than or equal to zero!".format(self.simSteps()))
                 break
         
             # Step forward in time
-            self.step(dt)
+            self.step(local_dt)
 
             #Print info
-            print_string = printer.getPrintString(t_end - self.simTime())
+            print_string = printer.getPrintString(self.simTime() - t_start)
             if (print_string):
                 self.logger.info("%s: %s", self, print_string)
                 try:

GPUSimulators/cuda/SWE2D_KP07.cu

@@ -201,12 +201,12 @@ __global__ void KP07Kernel(
         const int i = tx + 2; //Skip local ghost cells, i.e., +2
         const int j = ty + 2;
         
-        const float h1  = Q[0][j][i] + (F[0][ty][tx] - F[0][ty  ][tx+1]) * dt_ / dx_ 
-                                     + (G[0][ty][tx] - G[0][ty+1][tx  ]) * dt_ / dy_;
-        const float hu1 = Q[1][j][i] + (F[1][ty][tx] - F[1][ty  ][tx+1]) * dt_ / dx_ 
-                                     + (G[1][ty][tx] - G[1][ty+1][tx  ]) * dt_ / dy_;
-        const float hv1 = Q[2][j][i] + (F[2][ty][tx] - F[2][ty  ][tx+1]) * dt_ / dx_ 
-                                     + (G[2][ty][tx] - G[2][ty+1][tx  ]) * dt_ / dy_;
+        Q[0][j][i] += (F[0][ty][tx] - F[0][ty  ][tx+1]) * dt_ / dx_ 
+                    + (G[0][ty][tx] - G[0][ty+1][tx  ]) * dt_ / dy_;
+        Q[1][j][i] += (F[1][ty][tx] - F[1][ty  ][tx+1]) * dt_ / dx_ 
+                    + (G[1][ty][tx] - G[1][ty+1][tx  ]) * dt_ / dy_;
+        Q[2][j][i] += (F[2][ty][tx] - F[2][ty  ][tx+1]) * dt_ / dx_ 
+                    + (G[2][ty][tx] - G[2][ty+1][tx  ]) * dt_ / dy_;
 
         float* const h_row  = (float*) ((char*) h1_ptr_ + h1_pitch_*tj);
         float* const hu_row = (float*) ((char*) hu1_ptr_ + hu1_pitch_*tj);
@@ -214,14 +214,14 @@ __global__ void KP07Kernel(
 
         if (getOrder(step_order_) == 2 && getStep(step_order_) == 1) {
             //Write to main memory
-            h_row[ti]  = 0.5f*(h_row[ti]  + h1);
-            hu_row[ti] = 0.5f*(hu_row[ti] + hu1);
-            hv_row[ti] = 0.5f*(hv_row[ti] + hv1);
+            h_row[ti]  = 0.5f*(h_row[ti]  + Q[0][j][i]);
+            hu_row[ti] = 0.5f*(hu_row[ti] + Q[1][j][i]);
+            hv_row[ti] = 0.5f*(hv_row[ti] + Q[2][j][i]);
         }
         else {
-            h_row[ti] = h1;
-            hu_row[ti] = hu1;
-            hv_row[ti] = hv1;
+            h_row[ti]  = Q[0][j][i];
+            hu_row[ti] = Q[1][j][i];
+            hv_row[ti] = Q[2][j][i];
         }
     }
     

GPUSimulators/helpers/InitialConditions.py

@@ -22,6 +22,83 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 from GPUSimulators.Simulator import BoundaryCondition
 import numpy as np
+import gc
+
+
+
+def downsample(highres_solution, x_factor, y_factor=None):
+    if (y_factor == None):
+        y_factor = x_factor
+
+    assert(highres_solution.shape[1] % x_factor == 0)
+    assert(highres_solution.shape[0] % y_factor == 0)
+    
+    if (x_factor*y_factor == 1):
+        return highres_solution
+    
+    if (len(highres_solution.shape) == 1):
+        highres_solution = highres_solution.reshape((1, highres_solution.size))
+
+    nx = highres_solution.shape[1] / x_factor
+    ny = highres_solution.shape[0] / y_factor
+
+    return highres_solution.reshape([int(ny), int(y_factor), int(nx), int(x_factor)]).mean(3).mean(1)
+
+
+
+
+    
+def bump(nx, ny, width, height, bump_size=None, h_ref=0.5, h_amp=0.1, u_ref=0.0, u_amp=0.1, v_ref=0.0, v_amp=0.1, ref_nx=None, ref_ny=None):
+    
+    if (ref_nx == None):
+        ref_nx = nx
+    assert(ref_nx >= nx)
+      
+    if (ref_ny == None):
+        ref_ny = ny
+    assert(ref_ny >= ny)
+        
+    if (bump_size == None):
+        bump_size = width/5.0
+    
+    ref_dx = width / float(ref_nx)
+    ref_dy = height / float(ref_ny)
+
+    x_center = ref_dx*ref_nx/2.0
+    y_center = ref_dy*ref_ny/2.0
+    
+    x = ref_dx*(np.arange(0, ref_nx, dtype=np.float32)+0.5) - x_center
+    y = ref_dy*(np.arange(0, ref_ny, dtype=np.float32)+0.5) - y_center
+    xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
+    r = np.sqrt(xv**2 + yv**2)
+    xv = None
+    yv = None
+    gc.collect()
+    
+    #Generate highres then downsample
+    #h_highres = 0.5 + 0.1*np.exp(-(xv**2/size + yv**2/size))
+    h_highres = h_ref + h_amp*0.5*(1.0 + np.cos(np.pi*r/bump_size)) * (r < bump_size)
+    h = downsample(h_highres, ref_nx/nx, ref_ny/ny)
+    h_highres = None
+    gc.collect()
+    
+    #hu_highres = 0.1*np.exp(-(xv**2/size + yv**2/size))
+    u_highres = u_ref + u_amp*0.5*(1.0 + np.cos(np.pi*r/bump_size)) * (r < bump_size)
+    hu = downsample(u_highres, ref_nx/nx, ref_ny/ny)*h
+    u_highres = None
+    gc.collect()
+    
+    #hu_highres = 0.1*np.exp(-(xv**2/size + yv**2/size))
+    v_highres = v_ref + v_amp*0.5*(1.0 + np.cos(np.pi*r/bump_size)) * (r < bump_size)
+    hv = downsample(v_highres, ref_nx/nx, ref_ny/ny)*h
+    v_highres = None
+    gc.collect()
+    
+    dx = width/nx
+    dy = height/ny
+    
+    return h, hu, hv, dx, dy
+
 
 def genShockBubble(nx, ny, gamma):
     """
@@ -61,13 +138,8 @@ def genShockBubble(nx, ny, gamma):
     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,
@@ -80,7 +152,7 @@ def genShockBubble(nx, ny, gamma):
     arguments = {
         'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
         'nx': nx, 'ny': ny,
-        'dx': dx, 'dy': dy, 'dt': dt,
+        'dx': dx, 'dy': dy, 
         'g': g,
         'gamma': gamma,
         'boundary_conditions': bc
@@ -169,13 +241,8 @@ def genKelvinHelmholtz(nx, ny, gamma, roughness=0.125):
     
     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({
@@ -189,7 +256,7 @@ def genKelvinHelmholtz(nx, ny, gamma, roughness=0.125):
     arguments = {
         'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
         'nx': nx, 'ny': ny,
-        'dx': dx, 'dy': dy, 'dt': dt,
+        'dx': dx, 'dy': dy, 
         'g': g,
         'gamma': gamma,
         'boundary_conditions': bc
@@ -233,13 +300,8 @@ def genRayleighTaylor(nx, ny, gamma, version=0):
     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,
@@ -252,7 +314,7 @@ def genRayleighTaylor(nx, ny, gamma, version=0):
     arguments = {
         'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
         'nx': nx, 'ny': ny,
-        'dx': dx, 'dy': dy, 'dt': dt,
+        'dx': dx, 'dy': dy, 
         'g': g,
         'gamma': gamma,
         'boundary_conditions': bc