Merge pull request #471 from PrincetonUniversity/script_updates

Updated Python scripts

vis/python/athena_read.py

@@ -426,7 +426,7 @@ def athdf(filename, raw=False, data=None, quantities=None, dtype=None, level=Non
                             possible_max = (loc_this_dim+1) * 2**(level-level_this_dim)
                             num_blocks_this_dim = max(num_blocks_this_dim, possible_max)
                         else:
-                            possible_max = (loc_this_dim+1) / 2**(level_this_dim-level)
+                            possible_max = (loc_this_dim+1) // 2**(level_this_dim-level)
                             num_blocks_this_dim = max(num_blocks_this_dim, possible_max)
                     nx_vals.append(num_blocks_this_dim)
             elif block_size[d] == 1:  # singleton dimension
@@ -436,9 +436,9 @@ def athdf(filename, raw=False, data=None, quantities=None, dtype=None, level=Non
         nx1 = nx_vals[0]
         nx2 = nx_vals[1]
         nx3 = nx_vals[2]
-        lx1 = nx1 / block_size[0]
-        lx2 = nx2 / block_size[1]
-        lx3 = nx3 / block_size[2]
+        lx1 = nx1 // block_size[0]
+        lx2 = nx2 // block_size[1]
+        lx3 = nx3 // block_size[2]
         num_extended_dims = 0
         for nx in nx_vals:
             if nx > 1:
@@ -635,7 +635,7 @@ def center_func_3(xm, xp):
                     if np.all(levels == level):
                         data[xf] = np.empty(nx + 1, dtype=dtype)
                         for n_block in range(int((nx - 2*num_ghost)
-                                                 / (block_size[d-1] - 2*num_ghost))):
+                                                 // (block_size[d-1] - 2*num_ghost))):
                             sample_block = np.where(logical_locations[:, d-1]
                                                     == n_block)[0][0]
                             index_low = n_block * (block_size[d-1] - 2*num_ghost)
@@ -787,12 +787,12 @@ def center_func_3(xm, xp):
                 s = 2 ** (block_level - level)
 
                 # Calculate destination indices, without selection
-                il_d = block_location[0] * block_size[0] / s if nx1 > 1 else 0
-                jl_d = block_location[1] * block_size[1] / s if nx2 > 1 else 0
-                kl_d = block_location[2] * block_size[2] / s if nx3 > 1 else 0
-                iu_d = il_d + block_size[0] / s if nx1 > 1 else 1
-                ju_d = jl_d + block_size[1] / s if nx2 > 1 else 1
-                ku_d = kl_d + block_size[2] / s if nx3 > 1 else 1
+                il_d = block_location[0] * block_size[0] // s if nx1 > 1 else 0
+                jl_d = block_location[1] * block_size[1] // s if nx2 > 1 else 0
+                kl_d = block_location[2] * block_size[2] // s if nx3 > 1 else 0
+                iu_d = il_d + block_size[0] // s if nx1 > 1 else 1
+                ju_d = jl_d + block_size[1] // s if nx2 > 1 else 1
+                ku_d = kl_d + block_size[2] // s if nx3 > 1 else 1
 
                 # Calculate (restricted) source indices, with selection
                 il_s = max(il_d, i_min) - il_d
@@ -826,9 +826,9 @@ def center_func_3(xm, xp):
                 # Apply subsampling
                 if subsample:
                     # Calculate fine-level offsets (nearest cell at or below center)
-                    o1 = s/2 - 1 if nx1 > 1 else 0
-                    o2 = s/2 - 1 if nx2 > 1 else 0
-                    o3 = s/2 - 1 if nx3 > 1 else 0
+                    o1 = s//2 - 1 if nx1 > 1 else 0
+                    o2 = s//2 - 1 if nx2 > 1 else 0
+                    o3 = s//2 - 1 if nx3 > 1 else 0
 
                     # Assign values
                     for q, dataset, index in zip(quantities, quantity_datasets,
@@ -896,9 +896,9 @@ def center_func_3(xm, xp):
                                                                                block_num,
                                                                                k_s, j_s,
                                                                                i_s]
-                    loc1 = (nx1 > 1) * block_location[0] / s
-                    loc2 = (nx2 > 1) * block_location[1] / s
-                    loc3 = (nx3 > 1) * block_location[2] / s
+                    loc1 = (nx1 > 1) * block_location[0] // s
+                    loc2 = (nx2 > 1) * block_location[1] // s
+                    loc3 = (nx3 > 1) * block_location[2] // s
                     restricted_data[loc3, loc2, loc1] = True
 
             # Set level information for cells in this block
@@ -976,26 +976,27 @@ def restrict_like(vals, levels, vols=None):
         level_difference = max_level - level
         stride = 2 ** level_difference
         if nx3 > 1:
-            vals_level = np.reshape(vals * vols, (nx3/stride, stride, nx2/stride, stride,
-                                                  nx1/stride, stride))
-            vols_level = np.reshape(vols, (nx3/stride, stride, nx2/stride, stride,
-                                           nx1/stride, stride))
+            vals_level = np.reshape(vals * vols, (nx3//stride, stride, nx2//stride,
+                                                  stride, nx1//stride, stride))
+            vols_level = np.reshape(vols, (nx3//stride, stride, nx2//stride, stride,
+                                           nx1//stride, stride))
             vals_sum = np.sum(np.sum(np.sum(vals_level, axis=5), axis=3), axis=1)
             vols_sum = np.sum(np.sum(np.sum(vols_level, axis=5), axis=3), axis=1)
             vals_level = np.repeat(np.repeat(np.repeat(vals_sum / vols_sum, stride,
                                                        axis=0),
                                              stride, axis=1),
                                    stride, axis=2)
         elif nx2 > 1:
-            vals_level = np.reshape(vals * vols, (nx2/stride, stride, nx1/stride, stride))
-            vols_level = np.reshape(vols, (nx2/stride, stride, nx1/stride, stride))
+            vals_level = np.reshape(vals * vols, (nx2//stride, stride, nx1//stride,
+                                                  stride))
+            vols_level = np.reshape(vols, (nx2//stride, stride, nx1//stride, stride))
             vals_sum = np.sum(np.sum(vals_level, axis=3), axis=1)
             vols_sum = np.sum(np.sum(vols_level, axis=3), axis=1)
             vals_level = np.repeat(np.repeat(vals_sum / vols_sum, stride, axis=0),
                                    stride, axis=1)
         else:
-            vals_level = np.reshape(vals * vols, (nx1/stride, stride))
-            vols_level = np.reshape(vols, (nx1/stride, stride))
+            vals_level = np.reshape(vals * vols, (nx1//stride, stride))
+            vols_level = np.reshape(vols, (nx1//stride, stride))
             vals_sum = np.sum(vals_level, axis=1)
             vols_sum = np.sum(vols_level, axis=1)
             vals_level = np.repeat(vals_sum / vols_sum, stride, axis=0)

vis/python/plot_slice.py

@@ -75,9 +75,11 @@ def main(**kwargs):
     # Read data
     if quantities[0] == 'Levels':
         data = athena_read.athdf(kwargs['data_file'], quantities=quantities[1:],
-                                 level=level, return_levels=True, num_ghost=kwargs['num_ghost'])
+                                 level=level, return_levels=True,
+                                 num_ghost=kwargs['num_ghost'])
     else:
-        data = athena_read.athdf(kwargs['data_file'], quantities=quantities, level=level, num_ghost=kwargs['num_ghost'])
+        data = athena_read.athdf(kwargs['data_file'], quantities=quantities, level=level,
+                                 num_ghost=kwargs['num_ghost'])
 
     # Check that coordinates work with user choices
     coordinates = data['Coordinates'].decode('ascii', 'replace')

vis/python/plot_spherical.py

@@ -165,9 +165,9 @@ def theta_func(xmin, xmax, _, nf):
     # Perform slicing/averaging of scalar data
     if kwargs['midplane']:
         if nx2 % 2 == 0:
-            vals = np.mean(data[kwargs['quantity']][:, nx2/2-1:nx2/2+1, :], axis=1)
+            vals = np.mean(data[kwargs['quantity']][:, nx2//2-1:nx2//2+1, :], axis=1)
         else:
-            vals = data[kwargs['quantity']][:, nx2/2, :]
+            vals = data[kwargs['quantity']][:, nx2//2, :]
         if kwargs['average']:
             vals = np.repeat(np.mean(vals, axis=0, keepdims=True), nx3, axis=0)
     else:
@@ -177,8 +177,8 @@ def theta_func(xmin, xmax, _, nf):
         else:
             vals_right = 0.5 * (data[kwargs['quantity']]
                                 [-1, :, :] + data[kwargs['quantity']][0, :, :])
-            vals_left = 0.5 * (data[kwargs['quantity']][(nx3/2)-1, :, :]
-                               + data[kwargs['quantity']][nx3 / 2, :, :])
+            vals_left = 0.5 * (data[kwargs['quantity']][(nx3//2)-1, :, :]
+                               + data[kwargs['quantity']][nx3//2, :, :])
 
     # Join scalar data through boundaries
     if not kwargs['midplane']:
@@ -189,12 +189,12 @@ def theta_func(xmin, xmax, _, nf):
         if kwargs['midplane']:
             if nx2 % 2 == 0:
                 vals_r = np.mean(data[kwargs['stream'] + '1']
-                                 [:, nx2/2-1:nx2/2+1, :], axis=1).T
+                                 [:, nx2//2-1:nx2//2+1, :], axis=1).T
                 vals_phi = np.mean(data[kwargs['stream'] + '3']
-                                   [:, nx2/2-1:nx2/2+1, :], axis=1).T
+                                   [:, nx2//2-1:nx2//2+1, :], axis=1).T
             else:
-                vals_r = data[kwargs['stream'] + '1'][:, nx2/2, :].T
-                vals_phi = data[kwargs['stream'] + '3'][:, nx2/2, :].T
+                vals_r = data[kwargs['stream'] + '1'][:, nx2//2, :].T
+                vals_phi = data[kwargs['stream'] + '3'][:, nx2//2, :].T
             if kwargs['stream_average']:
                 vals_r = np.tile(np.reshape(np.mean(vals_r, axis=1), (nx1, 1)), nx3)
                 vals_phi = np.tile(np.reshape(np.mean(vals_phi, axis=1), (nx1, 1)), nx3)
@@ -206,9 +206,9 @@ def theta_func(xmin, xmax, _, nf):
                 vals_theta_left = -vals_theta_right
             else:
                 vals_r_right = data[kwargs['stream'] + '1'][0, :, :].T
-                vals_r_left = data[kwargs['stream'] + '1'][nx3/2, :, :].T
+                vals_r_left = data[kwargs['stream'] + '1'][nx3//2, :, :].T
                 vals_theta_right = data[kwargs['stream'] + '2'][0, :, :].T
-                vals_theta_left = -data[kwargs['stream'] + '2'][nx3/2, :, :].T
+                vals_theta_left = -data[kwargs['stream'] + '2'][nx3//2, :, :].T
 
     # Join vector data through boundaries
     if kwargs['stream'] is not None:
@@ -279,13 +279,13 @@ def theta_func(xmin, xmax, _, nf):
     vmin = kwargs['vmin']
     vmax = kwargs['vmax']
     if kwargs['logc']:
-        norm = colors.LogNorm()
+        norm = colors.LogNorm(vmin=vmin, vmax=vmax)
     else:
-        norm = colors.Normalize()
+        norm = colors.Normalize(vmin=vmin, vmax=vmax)
 
     # Make plot
     plt.figure()
-    im = plt.pcolormesh(x_grid, y_grid, vals, cmap=cmap, vmin=vmin, vmax=vmax, norm=norm)
+    im = plt.pcolormesh(x_grid, y_grid, vals, cmap=cmap, norm=norm)
     if kwargs['stream'] is not None:
         with warnings.catch_warnings():
             warnings.filterwarnings(

vis/python/spherical_refinement.py

@@ -115,8 +115,8 @@ def main(**kwargs):
     for L in range(max_levels+1):
 
         # Prepare description of all blocks at current level
-        num_blocks_r = num_r/num_r_block * 2**L
-        num_blocks_theta = num_theta/num_theta_block * 2**L
+        num_blocks_r = num_r//num_r_block * 2**L
+        num_blocks_theta = num_theta//num_theta_block * 2**L
         r_bounds.append(np.empty(num_blocks_r+1))
         theta_bounds.append(np.empty(num_blocks_theta+1))
         for i in range(num_blocks_r+1):
@@ -131,8 +131,8 @@ def main(**kwargs):
         refinement.append(np.ones((num_blocks_r, num_blocks_theta), dtype=bool))
         if L == 0:
             continue
-        for i in range(num_blocks_r/2):
-            for j in range(num_blocks_theta/2):
+        for i in range(num_blocks_r//2):
+            for j in range(num_blocks_theta//2):
                 if not refinement[L-1][i, j]:
                     refinement[L][i*2:(i+1)*2, j*2:(j+1)*2] = False
 
@@ -144,7 +144,7 @@ def main(**kwargs):
             for j in range(num_blocks_theta):
                 if not refinement[L][i, j]:
                     continue
-                if j < (num_blocks_theta+1)/2:
+                if j < (num_blocks_theta+1)//2:
                     theta1 = theta_adjust(theta_bounds[L][j], theta_compress)
                     theta2_unadjusted = pos_face(
                         theta_bounds[L][j], theta_bounds[L][j+1], 1.0, num_theta_block, 1)
@@ -158,11 +158,11 @@ def main(**kwargs):
                                              metric, parameters)
                 if min(w_r, w_theta, w_phi) < width_min:
                     refinement[L][i, j] = False
-                    refinement[L-1][i/2, j/2] = False
+                    refinement[L-1][i//2, j//2] = False
 
         # Make sure only entire blocks from previous level are refined
-        for i in range(num_blocks_r/2):
-            for j in range(num_blocks_theta/2):
+        for i in range(num_blocks_r//2):
+            for j in range(num_blocks_theta//2):
                 if not refinement[L-1][i, j]:
                     refinement[L][i*2:(i+1)*2, j*2:(j+1)*2] = False
 
@@ -185,22 +185,22 @@ def main(**kwargs):
     refinement_regions = []
     num_blocks = np.zeros(max_levels+1, dtype=int)
     num_blocks[0] = (
-        (num_r/num_r_block) * (num_theta/num_theta_block) * (num_phi/num_phi_block))
+        (num_r//num_r_block) * (num_theta//num_theta_block) * (num_phi//num_phi_block))
     for L in range(max_levels):
 
         # Calculate number of blocks in this level
-        num_blocks_r = num_r/num_r_block * 2**L
-        num_blocks_theta = num_theta/num_theta_block * 2**L
-        num_blocks_phi = num_phi/num_phi_block * 2**L
+        num_blocks_r = num_r//num_r_block * 2**L
+        num_blocks_theta = num_theta//num_theta_block * 2**L
+        num_blocks_phi = num_phi//num_phi_block * 2**L
 
         # Find block limit of region to be refined
         j_lims = np.empty(num_blocks_r, dtype=int)
         for i in range(num_blocks_r):
-            if not refinement[L][i, (num_blocks_theta-1)/2]:
+            if not refinement[L][i, (num_blocks_theta-1)//2]:
                 j_lims[i] = -1
                 continue
             j_lim = 0
-            for j in range((num_blocks_theta-1)/2, -1, -1):
+            for j in range((num_blocks_theta-1)//2, -1, -1):
                 if not refinement[L][i, j]:
                     j_lim = j+1
                     break
@@ -237,23 +237,23 @@ def main(**kwargs):
         print('\nThe following {} input blocks can be used to specify maximal refinement:'
               .format(num_refinement))
     for refinement_num, (L, i_start, i_end, j_lim) in enumerate(refinement_regions):
-        num_blocks_r = num_r/num_r_block * 2**L
-        num_blocks_theta = num_theta/num_theta_block * 2**L
+        num_blocks_r = num_r//num_r_block * 2**L
+        num_blocks_theta = num_theta//num_theta_block * 2**L
         if i_start == 0:
             r1 = r_min
         else:
             r1 = pos_face(r_bounds[L][i_start], r_bounds[L][i_start+1],
-                          r_ratio**(1.0/2**L), num_r_block, num_r_block/2)
+                          r_ratio**(1.0/2**L), num_r_block, num_r_block//2)
         if i_end == num_blocks_r-1:
             r2 = r_max
         else:
             r2 = pos_face(r_bounds[L][i_end], r_bounds[L][i_end+1], r_ratio**(1.0/2**L),
-                          num_r_block, num_r_block/2)
+                          num_r_block, num_r_block//2)
         if j_lim == 0:
             theta1 = theta_adjust(theta_min, theta_compress)
         else:
             theta1_unadjusted = pos_face(theta_bounds[L][j_lim], theta_bounds[L][j_lim+1],
-                                         1.0, num_theta_block, num_theta_block/2)
+                                         1.0, num_theta_block, num_theta_block//2)
             theta1 = theta_adjust(theta1_unadjusted, theta_compress)
         theta2 = np.pi - theta1
         print('\n<refinement{0}>'.format(refinement_num+1))

vis/python/uniform.py

@@ -28,7 +28,7 @@ def main(**kwargs):
         size = int(os.environ['OMPI_COMM_WORLD_SIZE'])
         rank = int(os.environ['OMPI_COMM_WORLD_RANK'])
         num_files = len(file_nums)
-        num_files_per_rank = num_files/size
+        num_files_per_rank = num_files//size
         num_files_extra = num_files % size
         num_files_list = ([num_files_per_rank+1] * num_files_extra
                           + [num_files_per_rank] * (size-num_files_extra))