Add UOST rotation

compass/ocean/tests/global_ocean/wave_mesh/__init__.py

@@ -67,7 +67,8 @@ def __init__(self, test_group, ocean_mesh, ocean_init):
         self.add_step(rotate_mesh_step)
 
         uost_file_step = WavesUostFiles(test_case=self,
-                                        wave_culled_mesh=cull_mesh_step)
+                                        wave_culled_mesh=cull_mesh_step,
+                                        wave_rotate_mesh=rotate_mesh_step)
         self.add_step(uost_file_step)
 
     def configure(self, config=None):

compass/ocean/tests/global_ocean/wave_mesh/uost_files.py

@@ -1,3 +1,4 @@
+import matplotlib.pyplot as plt
 import numpy as np
 
 from alphaBetaLab.alphaBetaLab.abEstimateAndSave import (
@@ -13,7 +14,7 @@ class WavesUostFiles(Step):
     """
     A step for creating the unresolved obstacles file for wave mesh
     """
-    def __init__(self, test_case, wave_culled_mesh,
+    def __init__(self, test_case, wave_culled_mesh, wave_rotate_mesh,
                  name='uost_files', subdir=None):
 
         super().__init__(test_case=test_case, name=name, subdir=subdir)
@@ -24,6 +25,11 @@ def __init__(self, test_case, wave_culled_mesh,
             filename='wave_mesh_culled.msh',
             work_dir_target=f'{culled_mesh_path}/wave_mesh_culled.msh')
 
+        angled_path = wave_rotate_mesh.path
+        self.add_input_file(
+            filename='angled.d',
+            work_dir_target=f'{angled_path}/angled.d')
+
         self.add_input_file(
             filename='etopo1_180.nc',
             target='etopo1_180.nc',
@@ -34,8 +40,10 @@ def run(self):
         Create unresolved obstacles for wave mesh and spectral resolution
         """
 
-        dirs = np.linspace(0, 2 * np.pi, 36)
-        nfreq = 50   # ET NOTE: this should be flexible
+        ndir = self.config.getint('wave_mesh', 'ndir')
+        nfreq = self.config.getint('wave_mesh', 'nfreq')
+
+        dirs = np.linspace(0, 2 * np.pi, ndir)
         minfrq = .035
         if (nfreq == 50):
             frqfactor = 1.07
@@ -50,20 +58,19 @@ def run(self):
         freqs = [minfrq * (frqfactor ** i) for i in range(1, nfreq + 1)]
 
         # definition of the spatial mesh
-        gridname = 'glo_unst'  # SB NOTE: should be flexible
+        gridname = 'waves_mesh_culled'  # SB NOTE: should be flexible
         mshfile = 'wave_mesh_culled.msh'
         triMeshSpec = triMeshSpecFromMshFile(mshfile)
 
         # path of the etopo1 bathymetry
-        # etopoFilePath = '/users/sbrus/scratch4/ \
-        # WW3_unstructured/GEBCO_2019.nc'
         etopoFilePath = 'etopo1_180.nc'
 
         # output directory
         outputDestDir = 'output/'
 
         # number of cores for parallel computing
-        nParWorker = 1  # SB NOTE: we should set this to the of cores on a node
+        nParWorker = 1
+        #nParWorker = self.config.getint('parallel', 'cores_per_node')
 
         # this option indicates that the computation
         # should be skipped for cells smaller than 3 km
@@ -77,3 +84,188 @@ def run(self):
         abEstimateAndSaveTriangularEtopo1(
             dirs, freqs, gridname, triMeshSpec, etopoFilePath,
             outputDestDir, nParWorker, abOptions=opt)
+
+        theta = np.radians(np.linspace(0.0, 360.0, ndir, endpoint=False))
+        freq = np.linspace(0.0, 1.0, nfreq)
+        Theta, Freq = np.meshgrid(theta, freq)
+
+        data = np.loadtxt('angled.d')
+        angled = data[:, 2]
+
+        filename_local_in = 'obstructions_local.waves_mesh_culled.in'
+        filename_local_out = 'obstructions_local.waves_mesh_culled.rtd.in'
+        lon_local, lat_local, nodes, sizes, alpha_local_avg, beta_local_avg, alpha_spec, beta_spec = self.read_alpha_beta(filename_local_in, nfreq)
+        alpha_interp, beta_interp = self.rotate_and_interpolate(Theta, nodes, angled, alpha_spec, beta_spec)
+        header = '$WAVEWATCH III LOCAL OBSTRUCTIONS'
+        self.write_alpha_beta(filename_local_out, header, nodes, lon_local, lat_local, sizes, alpha_interp, beta_interp)
+        self.plot_alpha_beta_spectra(Theta, Freq, alpha_spec, beta_spec, alpha_interp, beta_interp, angled, nodes, 'local')
+
+        filename_shadow_in = 'obstructions_shadow.waves_mesh_culled.in'
+        filename_shadow_out = 'obstructions_shadow.waves_mesh_culled.rtd.in'
+        lon_shadow, lat_shadow, nodes, sizes, alpha_shadow_avg, beta_shadow_avg, alpha_spec, beta_spec = self.read_alpha_beta(filename_shadow_in, nfreq)
+        alpha_interp, beta_interp = self.rotate_and_interpolate(Theta, nodes, angled, alpha_spec, beta_spec)
+        header = '$WAVEWATCH III SHADOW OBSTRUCTIONS'
+        self.write_alpha_beta(filename_shadow_out, header, nodes, lon_shadow, lat_shadow, sizes, alpha_interp, beta_interp)
+        self.plot_alpha_beta_spectra(Theta, Freq, alpha_spec, beta_spec, alpha_interp, beta_interp, angled, nodes, 'shadow')
+
+    def write_alpha_beta(self, filename, header, nodes, lon, lat, sizes, alpha_spec, beta_spec):
+
+        n = alpha_spec.shape[0]
+        nfreq = alpha_spec.shape[1]
+        ndir = alpha_spec.shape[2]
+
+        lines = []
+        lines.append(header)
+        lines.append(str(n))
+
+        for i in range(n):
+            lines.append('$ ilon ilat of the cell. lon: {:.8f}, lat: {:.8f}'.format(lon[i], lat[i]))
+            lines.append(str(nodes[i]) + '   1')
+            lines.append(sizes[i])
+            lines.append('$ mean alpha: {:.16}'.format(np.mean(alpha_spec[i, :, :])))
+            lines.append('$ mean beta: {:.16}'.format(np.mean(beta_spec[i, :, :])))
+            lines.append('$alpha by ik, ith')
+            for j in range(nfreq):
+                line = ''
+                for k in range(ndir):
+                    line = line + '{:.2f}  '.format(alpha_spec[i, j, k])
+                lines.append(line)
+            lines.append('$beta by ik, ith')
+            for j in range(nfreq):
+                line = ''
+                for k in range(ndir):
+                    line = line + '{:.2f}  '.format(beta_spec[i, j, k])
+                lines.append(line)
+
+        f = open(filename, 'w')
+        for line in lines:
+            f.write(line + '\n')
+        f.close()
+
+    def rotate_and_interpolate(self, Theta, nodes, angled,
+                               alpha_spec, beta_spec):
+
+        n = alpha_spec.shape[0]
+        nfreq = alpha_spec.shape[1]
+        ndir = alpha_spec.shape[2]
+
+        alpha_interp = np.zeros((n, nfreq, ndir))
+        beta_interp = np.zeros((n, nfreq, ndir))
+        for i in range(n):
+            nd = nodes[i] - 1
+            Theta2 = Theta - angled[nd] * np.pi / 180.0
+            for j in range(nfreq):
+                alpha_interp[i, j, :] = np.interp(Theta2[j, :],
+                                                  Theta[j, :],
+                                                  alpha_spec[i, j, :],
+                                                  period=2.0 * np.pi)
+                beta_interp[i, j, :] = np.interp(Theta2[j, :],
+                                                 Theta[j, :],
+                                                 beta_spec[i, j, :],
+                                                 period=2.0 * np.pi)
+
+        return [alpha_interp, beta_interp]
+
+    def plot_alpha_beta_spectra(self, Theta, Freq, alpha_spec, beta_spec,
+                                alpha_interp, beta_interp, angled, nodes, kind):
+
+        for i in range(10):
+            print(i)
+
+            fig = plt.figure(figsize=[8, 4])
+            ax = fig.add_subplot(2, 2, 1, polar=True)
+            ax.set_theta_zero_location("N")
+            ax.set_theta_direction(-1)
+            ax.contourf(Theta, Freq, alpha_spec[i, :, :], 30)
+
+            ax = fig.add_subplot(2, 2, 2, polar=True)
+            ax.set_theta_zero_location("N")
+            ax.set_theta_direction(-1)
+            ax.contourf(Theta, Freq, beta_spec[i, :, :], 30)
+
+            ax = fig.add_subplot(2, 2, 3, polar=True)
+            ax.set_theta_zero_location("N")
+            ax.set_theta_direction(-1)
+            ax.contourf(Theta, Freq, alpha_interp[i, :, :], 30)
+            ax.set_title('AnglD = ' + str(angled[nodes[i] - 1]))
+
+            ax = fig.add_subplot(2, 2, 4, polar=True)
+            ax.set_theta_zero_location("N")
+            ax.set_theta_direction(-1)
+            ax.contourf(Theta, Freq, beta_interp[i, :, :], 30)
+            ax.set_title('AnglD = ' + str(angled[nodes[i] - 1]))
+
+            plt.savefig(kind + '_spec_' + str(i) + '.png', bbox_inches='tight')
+
+    def read_alpha_beta(self, filename, nfreq):
+        f = open(filename, 'r')
+        lines = f.read().splitlines()
+
+        nodes = []
+        lon = []
+        lat = []
+        sizes = []
+        alpha_avg = []
+        beta_avg = []
+        alpha_spec = []
+        beta_spec = []
+
+        line = 1  # header comment
+        n = int(lines[line])
+        for i in range(n):
+
+            line = line + 1  # lon lat comment
+
+            text = lines[line]
+            text_sp = text.split()
+            x = float(text_sp[7].replace(',', ''))
+            y = float(text_sp[9])
+
+            line = line + 1  # node number
+            nodes.append(int(lines[line].split()[0]))
+
+            line = line + 1  # sizes comment
+            line = line + 1  # sizes
+            sizes.append(lines[line])
+
+            line = line + 1  # mean alpha
+            text = lines[line]
+            text_sp = text.split()
+            a = float(text_sp[-1])
+
+            line = line + 1  # mean beta
+            text = lines[line]
+            text_sp = text.split()
+            b = float(text_sp[-1])
+
+            line = line + 1  # alpha comment
+            spectrum = []
+            for i in range(nfreq):
+                line = line + 1
+                spectrum.append(lines[line].split())
+            alpha_spec.append(spectrum)
+            del spectrum
+
+            line = line + 1  # beta comment
+            spectrum = []
+            for i in range(nfreq):
+                line = line + 1
+                spectrum.append(lines[line].split())
+            beta_spec.append(spectrum)
+            del spectrum
+
+            lon.append(x)
+            lat.append(y)
+            alpha_avg.append(a)
+            beta_avg.append(b)
+
+        lon = np.array(lon)
+        lat = np.array(lat)
+        nodes = np.array(nodes)
+        alpha_avg = np.array(alpha_avg)
+        beta_avg = np.array(beta_avg)
+        alpha_spec = np.array(alpha_spec)
+        beta_spec = np.array(beta_spec)
+
+        return [lon, lat, nodes, sizes, alpha_avg, beta_avg,
+                alpha_spec, beta_spec]

compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg

@@ -8,3 +8,6 @@ depth_threshold_global = 1000.0
 distance_threshold_global = 300.0
 refined_res = 20000.0
 maxres = 225000.0
+
+ndir = 36
+nfreq = 50