From fdac082adc907ff04c6f52baba3b9c4d2f56e41c Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Mon, 15 Apr 2024 09:49:29 -0700 Subject: [PATCH 01/26] Start wave mesh test case --- compass/ocean/tests/global_ocean/__init__.py | 11 +++- .../tests/global_ocean/wave_mesh/__init__.py | 55 +++++++++++++++++++ .../tests/global_ocean/wave_mesh/base_mesh.py | 38 +++++++++++++ 3 files changed, 103 insertions(+), 1 deletion(-) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/__init__.py create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index 3d266a9cc9..de4244626d 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -1,3 +1,5 @@ +from comass.ocean.test.global_ocean.wave_mesh import WaveMesh + from compass.ocean.tests.global_ocean.analysis_test import AnalysisTest from compass.ocean.tests.global_ocean.daily_output_test import DailyOutputTest from compass.ocean.tests.global_ocean.decomp_test import DecompTest @@ -62,10 +64,14 @@ def __init__(self, mpas_core): # A test case for making E3SM support files from an existing mesh self.add_test_case(FilesForE3SM(test_group=self)) + # Create waves mesh + self._add_tests(mesh_names=['Icos'], include_wave_mesh=True) + def _add_tests(self, mesh_names, high_res_topography=True, include_rk4=False, include_regression=False, include_phc=False, - include_en4_1900=False): + include_en4_1900=False, + include_wave_mesh=False): """ Add test cases for the given mesh(es) """ default_ic = 'WOA23' @@ -167,3 +173,6 @@ def _add_tests(self, mesh_names, high_res_topography=True, FilesForE3SM( test_group=self, mesh=mesh_test, init=init_test, dynamic_adjustment=dynamic_adjustment_test)) + + if include_wave_mesh: + self.add_test_case(WaveMesh(test_group=self, ocean_mesh=mest_test)) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py new file mode 100644 index 0000000000..c1d91f1f77 --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -0,0 +1,55 @@ +#from compass.ocean.tests.global_ocean.metadata import ( +# get_author_and_email_from_git, +#) +from compass.testcase import TestCase + +#from compass.validate import compare_variables + + +class WaveMesh(TestCase): + """ + A test case for creating a global MPAS-Ocean mesh + + Attributes + ---------- + """ + def __init__(self, test_group, ocean_mesh): + """ + Create test case for creating a global MPAS-Ocean mesh + + Parameters + ---------- + test_group : compass.ocean.tests.global_ocean.GlobalOcean + The global ocean test group that this test case belongs to + + """ + name = 'wave_mesh' + subdir = f'{mesh_name}/{name}' + super().__init__(test_group=test_group, name=name, subdir=subdir) + + self.package = f'compass.ocean.tests.global_ocean.mesh.{mesh_lower}' + self.mesh_config_filename = f'{mesh_lower}.cfg' + + self.add_step(base_mesh_step) + +# def configure(self, config=None): +# """ +# Modify the configuration options for this test case +# +# config : compass.config.CompassConfigParser, optional +# Configuration options to update if not those for this test case +# """ +# if config is None: +# config = self.config +# config.add_from_package('compass.mesh', 'mesh.cfg', exception=True) +# +# get_author_and_email_from_git(config) + +# def validate(self): +# """ +# Test cases can override this method to perform validation of variables +# and timers +# """ +# variables = ['xCell', 'yCell', 'zCell'] +# compare_variables(test_case=self, variables=variables, +# filename1='cull_mesh/culled_mesh.nc') diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py new file mode 100644 index 0000000000..e294f6b45b --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -0,0 +1,38 @@ +import numpy as np + +from compass.mesh import QuasiUniformSphericalMeshStep + + +class WavesBaseMesh(QuasiUniformSphericalMeshStep): + """ + A step for creating wave mesh based on an ocean mesh + """ + pass + #def build_cell_width_lat_lon(self): + # """ + # Create cell width array for this mesh on a regular latitude-longitude + # grid + + # Returns + # ------- + # cellWidth : numpy.array + # m x n array of cell width in km + + # lon : numpy.array + # longitude in degrees (length n and between -180 and 180) + + # lat : numpy.array + # longitude in degrees (length m and between -90 and 90) + # """ + + # dlon = 10. + # dlat = 0.1 + # nlon = int(360. / dlon) + 1 + # nlat = int(180. / dlat) + 1 + # lon = np.linspace(-180., 180., nlon) + # lat = np.linspace(-90., 90., nlat) + + # cellWidthVsLat = mdt.EC_CellWidthVsLat(lat) + # cellWidth = np.outer(cellWidthVsLat, np.ones([1, lon.size])) + + # return cellWidth, lon, lat From e4a5e8793ca96f24033629a4c85863d3f6c61fb1 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Mon, 15 Apr 2024 13:43:00 -0500 Subject: [PATCH 02/26] Fix wave_mesh test case typos --- compass/ocean/tests/global_ocean/__init__.py | 6 +++--- .../ocean/tests/global_ocean/wave_mesh/__init__.py | 13 ++++++++----- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index de4244626d..f11572bc89 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -1,5 +1,3 @@ -from comass.ocean.test.global_ocean.wave_mesh import WaveMesh - from compass.ocean.tests.global_ocean.analysis_test import AnalysisTest from compass.ocean.tests.global_ocean.daily_output_test import DailyOutputTest from compass.ocean.tests.global_ocean.decomp_test import DecompTest @@ -15,6 +13,7 @@ from compass.ocean.tests.global_ocean.performance_test import PerformanceTest from compass.ocean.tests.global_ocean.restart_test import RestartTest from compass.ocean.tests.global_ocean.threads_test import ThreadsTest +from compass.ocean.tests.global_ocean.wave_mesh import WaveMesh from compass.testgroup import TestGroup @@ -175,4 +174,5 @@ def _add_tests(self, mesh_names, high_res_topography=True, dynamic_adjustment=dynamic_adjustment_test)) if include_wave_mesh: - self.add_test_case(WaveMesh(test_group=self, ocean_mesh=mest_test)) + self.add_test_case( + WaveMesh(test_group=self, ocean_mesh=mesh_test)) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index c1d91f1f77..7f6cf0e8d3 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -1,9 +1,10 @@ -#from compass.ocean.tests.global_ocean.metadata import ( -# get_author_and_email_from_git, -#) +# from compass.ocean.tests.global_ocean.metadata import ( +# get_author_and_email_from_git, +# ) +from compass.ocean.tests.global_ocean.wave_mesh.base_mesh import WavesBaseMesh from compass.testcase import TestCase -#from compass.validate import compare_variables +# from compass.validate import compare_variables class WaveMesh(TestCase): @@ -24,12 +25,14 @@ def __init__(self, test_group, ocean_mesh): """ name = 'wave_mesh' - subdir = f'{mesh_name}/{name}' + subdir = f'{ocean_mesh.mesh_name}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) + mesh_lower = ocean_mesh.mesh_name.lower() self.package = f'compass.ocean.tests.global_ocean.mesh.{mesh_lower}' self.mesh_config_filename = f'{mesh_lower}.cfg' + base_mesh_step = WavesBaseMesh(test_case=self) self.add_step(base_mesh_step) # def configure(self, config=None): From 08bde7584f8c279681f5adb98b2a730adf6c0ac3 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Wed, 17 Apr 2024 12:06:28 -0500 Subject: [PATCH 03/26] Update base_mesh step for waves --- .../tests/global_ocean/wave_mesh/__init__.py | 32 +- .../tests/global_ocean/wave_mesh/base_mesh.py | 366 ++++++++++++++++-- .../global_ocean/wave_mesh/wave_mesh.cfg | 10 + 3 files changed, 363 insertions(+), 45 deletions(-) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 7f6cf0e8d3..cb59bdb214 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -28,25 +28,25 @@ def __init__(self, test_group, ocean_mesh): subdir = f'{ocean_mesh.mesh_name}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) - mesh_lower = ocean_mesh.mesh_name.lower() - self.package = f'compass.ocean.tests.global_ocean.mesh.{mesh_lower}' - self.mesh_config_filename = f'{mesh_lower}.cfg' + self.package = 'compass.ocean.tests.global_ocean.wave_mesh' + self.mesh_config_filename = 'wave_mesh.cfg' - base_mesh_step = WavesBaseMesh(test_case=self) + base_mesh_step = WavesBaseMesh(test_case=self, + ocean_mesh=ocean_mesh) self.add_step(base_mesh_step) -# def configure(self, config=None): -# """ -# Modify the configuration options for this test case -# -# config : compass.config.CompassConfigParser, optional -# Configuration options to update if not those for this test case -# """ -# if config is None: -# config = self.config -# config.add_from_package('compass.mesh', 'mesh.cfg', exception=True) -# -# get_author_and_email_from_git(config) + def configure(self, config=None): + """ + Modify the configuration options for this test case + + config : compass.config.CompassConfigParser, optional + Configuration options to update if not those for this test case + """ + if config is None: + config = self.config + config.add_from_package('compass.mesh', 'mesh.cfg', exception=True) + + # get_author_and_email_from_git(config) # def validate(self): # """ diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index e294f6b45b..2bdc994b04 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -1,4 +1,15 @@ +import timeit + +import cartopy # noqa: F401 +import cartopy.crs as ccrs # noqa: F401 +import cartopy.feature as cfeature # noqa: F401 +import jigsawpy +import matplotlib.pyplot as plt import numpy as np +import shapefile +from mpas_tools.cime.constants import constants +from netCDF4 import Dataset +from scipy import interpolate, spatial from compass.mesh import QuasiUniformSphericalMeshStep @@ -7,32 +18,329 @@ class WavesBaseMesh(QuasiUniformSphericalMeshStep): """ A step for creating wave mesh based on an ocean mesh """ - pass - #def build_cell_width_lat_lon(self): - # """ - # Create cell width array for this mesh on a regular latitude-longitude - # grid - - # Returns - # ------- - # cellWidth : numpy.array - # m x n array of cell width in km - - # lon : numpy.array - # longitude in degrees (length n and between -180 and 180) - - # lat : numpy.array - # longitude in degrees (length m and between -90 and 90) - # """ - - # dlon = 10. - # dlat = 0.1 - # nlon = int(360. / dlon) + 1 - # nlat = int(180. / dlat) + 1 - # lon = np.linspace(-180., 180., nlon) - # lat = np.linspace(-90., 90., nlat) - - # cellWidthVsLat = mdt.EC_CellWidthVsLat(lat) - # cellWidth = np.outer(cellWidthVsLat, np.ones([1, lon.size])) - - # return cellWidth, lon, lat + def __init__(self, test_case, ocean_mesh, name='base_mesh', subdir=None): + + super().__init__(test_case=test_case, name=name, subdir=subdir, + cell_width=None) + + base_mesh_path = ocean_mesh.steps['base_mesh'].path + self.add_input_file( + filename='ocean_mesh.nc', + work_dir_target=f'{base_mesh_path}/base_mesh.nc') + + def setup(self): + + super().setup() + self.opts.init_file = 'init.msh' + + def build_cell_width_lat_lon(self): + """ + Create cell width array for this mesh on a regular latitude-longitude + grid + + Returns + ------- + cellWidth : numpy.array + m x n array of cell width in km + + lon : numpy.array + longitude in degrees (length n and between -180 and 180) + + lat : numpy.array + longitude in degrees (length m and between -90 and 90) + """ + km = 1000.0 + + lon_min = -180.0 + lon_max = 180.0 + dlon = self.config.getfloat('wave_mesh', 'hfun_grid_spacing') + nlon = int((lon_max - lon_min) / dlon) + 1 + lat_min = -90.0 + lat_max = 90.0 + nlat = int((lat_max - lat_min) / dlon) + 1 + + xlon = np.linspace(lon_min, lon_max, nlon) + ylat = np.linspace(lat_min, lat_max, nlat) + + earth_radius = constants['SHR_CONST_REARTH'] / km + shapefiles = ['GSHHS_l_L1.shp', 'GSHHS_l_L6.shp'] + cell_width = self.cell_widthVsLatLon(xlon, ylat, shapefiles, + earth_radius, 'ocean_mesh.nc') + cell_width = cell_width / km + + self.create_initial_points('ocean_mesh.nc', xlon, ylat, cell_width, + earth_radius, self.opts.init_file) + + return cell_width, xlon, ylat + + def cell_widthVsLatLon(self, lon, lat, coast_shapefiles, + sphere_radius, ocean_mesh): + + config = self.config + + depth_threshold_refined = config.getfloat('wave_mesh', + 'depth_threshold_refined') + dist_threshold_refined = config.getfloat('wave_mesh', + 'distance_threshold_refined') + depth_threshold_global = config.getfloat('wave_mesh', + 'depth_threshold_global') + dist_threshold_global = config.getfloat('wave_mesh', + 'distance_threshold_global') + refined_res = config.getfloat('wave_mesh', 'refined_res') + maxres = config.getfloat('wave_mesh', 'maxres') + + # Create structrued grid points + nlon = lon.size + nlat = lat.size + Lon_grd, Lat_grd = np.meshgrid(np.radians(lon), np.radians(lat)) + xy_pts = np.vstack((Lon_grd.ravel(), Lat_grd.ravel())).T + + # Create structured grid points and background cell_with array + cell_width = np.zeros(Lon_grd.shape) + maxres + + # Get ocean mesh variables + nc_file = Dataset(ocean_mesh, 'r') + areaCell = nc_file.variables['areaCell'][:] + lonCell = nc_file.variables['lonCell'][:] + latCell = nc_file.variables['latCell'][:] + bottomDepth = nc_file.variables['bottomDepth'][:] + + # Transform 0,360 range to -180,180 + idx, = np.where(lonCell > np.pi) + lonCell[idx] = lonCell[idx] - 2.0 * np.pi + idx, = np.where(lonCell < -np.pi) + lonCell[idx] = lonCell[idx] + 2.0 * np.pi + + # Interpolate cellWidth onto background grid + cellWidth = 2.0 * np.sqrt(areaCell / np.pi) + hfun = interpolate.NearestNDInterpolator((lonCell, latCell), cellWidth) + hfun_interp = hfun(xy_pts) + hfun_grd = np.reshape(hfun_interp, (nlat, nlon)) + + # Interpolate bathymetry onto background grid + bathy = interpolate.NearestNDInterpolator((lonCell, latCell), + bottomDepth) + bathy_interp = bathy(xy_pts) + bathy_grd = np.reshape(bathy_interp, (nlat, nlon)) + + # Get distance to coasts + D = self.distance_to_shapefile_points(coast_shapefiles, lon, lat, + sphere_radius, reggrid=True) + + # Apply refined region criteria + idxx, idxy = np.where((bathy_grd < depth_threshold_refined) & + (bathy_grd > 0.0) & + (D < dist_threshold_refined) & + (hfun_grd < refined_res)) + cell_width[idxx, idxy] = hfun_grd[idxx, idxy] + + # Apply global region criteria + idxx, idxy = np.where((bathy_grd < depth_threshold_global) & + (bathy_grd > 0.0) & + (D < dist_threshold_global)) + cell_width[idxx, idxy] = hfun_grd[idxx, idxy] + + hfun_slope_lim = config.getfloat('wave_mesh', 'hfun_slope_lim') + self.limit_spacing_gradient(lon, lat, cell_width, hfun_slope_lim) + + # Plot + fig = plt.figure(figsize=(16, 8)) + + ax1 = fig.add_subplot(4, 1, 1, projection=ccrs.PlateCarree()) + plt1 = ax1.contourf(lon, lat, cell_width, transform=ccrs.PlateCarree()) + ax1.set_title('waves mesh cell width') + fig.colorbar(plt1, ax=ax1) + + ax2 = fig.add_subplot(4, 1, 2, projection=ccrs.PlateCarree()) + plt2 = ax2.contourf(lon, lat, hfun_grd, transform=ccrs.PlateCarree()) + ax2.set_title('ocean cell width') + fig.colorbar(plt2, ax=ax2) + + ax3 = fig.add_subplot(4, 1, 3, projection=ccrs.PlateCarree()) + plt3 = ax3.contourf(lon, lat, bathy_grd, transform=ccrs.PlateCarree()) + ax3.set_title('bathymetry') + fig.colorbar(plt3, ax=ax3) + + ax4 = fig.add_subplot(4, 1, 4, projection=ccrs.PlateCarree()) + plt4 = ax4.contourf(lon, lat, D, transform=ccrs.PlateCarree()) + ax4.set_title('distance to coast') + fig.colorbar(plt4, ax=ax4) + + fig.savefig('cell_width.png', bbox_inches='tight', dpi=400) + plt.close() + + return cell_width + + def limit_spacing_gradient(self, lon, lat, cell_width, dhdx): + + print("Smoothing h(x) via |dh/dx| limits...") + + opts = jigsawpy.jigsaw_jig_t() + opts.hfun_file = "spac_pre_smooth.msh" + opts.jcfg_file = "opts_pre_smooth.jig" + opts.verbosity = +1 + + spac = jigsawpy.jigsaw_msh_t() + spac.mshID = "ellipsoid-grid" + spac.xgrid = lon + spac.ygrid = lat + spac.value = cell_width + spac.slope = np.full(spac.value.shape, dhdx) + jigsawpy.savemsh(opts.hfun_file, spac) + + jigsawpy.cmd.marche(opts, spac) + + return spac.value + + def distance_to_shapefile_points(self, shpfiles, lon, lat, + sphere_radius, reggrid=False): + + # Get coastline coordinates from shapefile + pt_list = [] + for shpfile in shpfiles: + + sf = shapefile.Reader(shpfile) + shapes = sf.shapes() + # records = sf.records() + n = len(sf.shapes()) + + for i in range(n): + + points = shapes[i].points + + if len(points) < 20: + continue + + for pt in points: + pt_list.append([pt[0], pt[1]]) + + coast_pts = np.radians(np.array(pt_list)) + + # Convert coastline points to x,y,z and create kd-tree + npts = coast_pts.shape[0] + coast_pts_xyz = np.zeros((npts, 3)) + x, y, z = self.lonlat2xyz(coast_pts[:, 0], coast_pts[:, 1], + sphere_radius) + coast_pts_xyz[:, 0] = x + coast_pts_xyz[:, 1] = y + coast_pts_xyz[:, 2] = z + tree = spatial.KDTree(coast_pts_xyz) + + # Make sure longitude and latitude are in radians + if np.amax(lon) < 2.1 * np.pi: + lonr = lon + latr = lat + else: + lonr = np.radians(lon) + latr = np.radians(lat) + + # Convert backgound grid coordinates to x,y,z + # and put in a nx x 3 array for kd-tree query + if reggrid: + Lon, Lat = np.meshgrid(lonr, latr) + else: + Lon = lonr + Lat = latr + X, Y, Z = self.lonlat2xyz(Lon, Lat, sphere_radius) + pts = np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T + + # Find distances of background grid coordinates to the coast + print(" Finding distance") + start = timeit.default_timer() + d, idx = tree.query(pts) + end = timeit.default_timer() + print(" Done") + print(" " + str(end - start) + " seconds") + + D = np.reshape(d, Lon.shape) + + return D + + def create_initial_points(self, meshfile, lon, lat, hfunction, + sphere_radius, outfile): + + # Open MPAS mesh and get cell variables + nc_file = Dataset(meshfile, 'r') + lonCell = nc_file.variables['lonCell'][:] + latCell = nc_file.variables['latCell'][:] + nCells = lonCell.shape[0] + + # Transform 0,360 range to -180,180 + idx, = np.where(lonCell > np.pi) + lonCell[idx] = lonCell[idx] - 2.0 * np.pi + + # Interpolate hfunction onto mesh cell centers + hfun = interpolate.RegularGridInterpolator( + (np.radians(lon), np.radians(lat)), + hfunction.T) + mesh_pts = np.vstack((lonCell, latCell)).T + hfun_interp = hfun(mesh_pts) + + # Find cells in refined region of waves mesh + max_res = np.amax(hfunction) + idx, = np.where(hfun_interp < 0.5 * max_res) + + # Find boundary cells + # some meshes have non-zero values + # in the extra columns of cellsOnCell + # in this case, these must be zeroed out + # to correctly identify boundary cells + nEdgesOnCell = nc_file.variables['nEdgesOnCell'][:] + cellsOnCell = nc_file.variables['cellsOnCell'][:] + nz = np.zeros(cellsOnCell.shape, dtype=bool) + for i in range(nCells): + nz[i, 0:nEdgesOnCell[i]] = True + cellsOnCell[~nz] = 0 + nCellsOnCell = np.count_nonzero(cellsOnCell, axis=1) + is_boundary_cell = np.equal(nCellsOnCell, nEdgesOnCell) + idx_bnd, = np.where(is_boundary_cell == False) # noqa: E712 + + # Force inclusion of all boundary cells + idx = np.union1d(idx, idx_bnd) + + # Get coordinates of points + lon = lonCell[idx] + lat = latCell[idx] + + lon = np.append(lon, 0.0) + lat = np.append(lat, 0.5 * np.pi) + + npt = lon.size + print(npt) + + # Change to Cartesian coordinates + x, y, z = self.lonlat2xyz(lon, lat, sphere_radius) + + # Get coordinates and ID into structured array + # (for use with np.savetxt) + pt_list = [] + for i in range(npt): + # ID of -1 specifies that node is fixed + pt_list.append((x[i], y[i], z[i], -1)) + pt_type = np.dtype({'names': ['x', 'y', 'z', 'id'], + 'formats': [np.float64, np.float64, + np.float64, np.int32]}) + pts = np.array(pt_list, dtype=pt_type) + + # Write initial conditions file + f = open(outfile, 'w') + f.write('# Initial coordinates \n') + f.write('MSHID=3;EUCLIDEAN-MESH\n') + f.write('NDIMS=3\n') + f.write(f'POINT={npt}\n') + np.savetxt(f, pts, fmt='%.12e;%.12e;%.12e;%2i') + f.close() + + init = jigsawpy.jigsaw_msh_t() + jigsawpy.loadmsh(self.opts.init_file, init) + jigsawpy.savevtk("init.vtk", init) + + return + + def lonlat2xyz(self, lon, lat, R): + + x = R * np.multiply(np.cos(lon), np.cos(lat)) + y = R * np.multiply(np.sin(lon), np.cos(lat)) + z = R * np.sin(lat) + + return x, y, z diff --git a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg new file mode 100644 index 0000000000..2ecf947586 --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg @@ -0,0 +1,10 @@ +[wave_mesh] + +hfun_grid_spacing = 0.5 +hfun_slope_lim = 0.15 +depth_threshold_refined = 1000.0 +distance_threshold_refined = 300.0 +depth_threshold_global = 1000.0 +distance_threshold_global = 300.0 +refined_res = 20000.0 +maxres = 225000.0 From d3fe8d5482c8436648a76ee3cada53b3d38525bd Mon Sep 17 00:00:00 2001 From: Erin Thomas Date: Wed, 17 Apr 2024 14:50:24 -0500 Subject: [PATCH 04/26] uost rebase --- .../tests/global_ocean/wave_mesh/__init__.py | 6 ++ .../global_ocean/wave_mesh/uost_files.py | 100 ++++++++++++++++++ 2 files changed, 106 insertions(+) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/uost_files.py diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index cb59bdb214..0d550c8cc4 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -2,6 +2,9 @@ # get_author_and_email_from_git, # ) from compass.ocean.tests.global_ocean.wave_mesh.base_mesh import WavesBaseMesh +from compass.ocean.tests.global_ocean.wave_mesh.uost_files import ( + WavesUostFiles, +) from compass.testcase import TestCase # from compass.validate import compare_variables @@ -34,6 +37,9 @@ def __init__(self, test_group, ocean_mesh): base_mesh_step = WavesBaseMesh(test_case=self, ocean_mesh=ocean_mesh) self.add_step(base_mesh_step) + + uost_file_step = WavesUostFiles(test_case=self) + self.add_step(uost_file_step) def configure(self, config=None): """ diff --git a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py new file mode 100644 index 0000000000..b5a758584b --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py @@ -0,0 +1,100 @@ +import numpy as np + +from alphaBetaLab.alphaBetaLab.abEstimateAndSave import ( + abEstimateAndSaveTriangularEtopo1, + triMeshSpecFromMshFile, +) +# importing from alphaBetaLab the needed components +from alphaBetaLab.alphaBetaLab.abOptionManager import abOptions + +# from compass.mesh import QuasiUniformSphericalMeshStep + + +class WavesUostFiles(): + """ + A step for creating the unresolved obstacles file for wave mesh + """ + # pass + + def __init__(self, test_case): + """ + Create a new step + """ + super().__init__(test_case, name='uost_files') + + def run(self): + """ + Run this step + """ + super().run() + + dirs = np.linspace(0, 2 * np.pi, 36) + nfreq = 50 # ET NOTE: this should be flexible + minfrq = .035 + if (nfreq == 50): + frqfactor = 1.07 + elif (nfreq == 36): + frqfactor = 1.10 + elif (nfreq == 25): + frqfactor = 1.147 + # PRINT SOME ERROR MESSAGE FOR NON-supoorted spectral reslution + + freqs = [minfrq * (frqfactor ** i) for i in range(1, nfreq + 1)] + + # definition of the spatial mesh + gridname = 'glo_unst' + mshfile = 'global.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 + + # this option indicates that the computation + # should be skipped for cells smaller than 3 km + minSizeKm = 3 + opt = abOptions(minSizeKm=minSizeKm) + + # instruction to do the computation and save the output + # abEstimateAndSaveTriangularGebco( + # dirs, freqs, gridname, triMeshSpec, + # etopoFilePath, outputDestDir, nParWorker, abOptions=opt) + abEstimateAndSaveTriangularEtopo1( + dirs, freqs, gridname, triMeshSpec, etopoFilePath, + outputDestDir, nParWorker, abOptions=opt) + + # def build_cell_width_lat_lon(self): + # """ + # Create cell width array for this mesh on a regular latitude-longitude + # grid + + # Returns + # ------- + # cellWidth : numpy.array + # m x n array of cell width in km + + # lon : numpy.array + # longitude in degrees (length n and between -180 and 180) + + # lat : numpy.array + # longitude in degrees (length m and between -90 and 90) + # """ + + # dlon = 10. + # dlat = 0.1 + # nlon = int(360. / dlon) + 1 + # nlat = int(180. / dlat) + 1 + # lon = np.linspace(-180., 180., nlon) + # lat = np.linspace(-90., 90., nlat) + + # cellWidthVsLat = mdt.EC_CellWidthVsLat(lat) + # cellWidth = np.outer(cellWidthVsLat, np.ones([1, lon.size])) + + # return cellWidth, lon, lat From 591e54d9d90176f5a7041e28c42be0d8c9f7ce24 Mon Sep 17 00:00:00 2001 From: Erin Thomas Date: Wed, 17 Apr 2024 15:40:42 -0500 Subject: [PATCH 05/26] UOST files step --- .../global_ocean/wave_mesh/uost_files.py | 61 ++++++------------- 1 file changed, 20 insertions(+), 41 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py index b5a758584b..e0e99bbacc 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py @@ -6,27 +6,33 @@ ) # importing from alphaBetaLab the needed components from alphaBetaLab.alphaBetaLab.abOptionManager import abOptions +from compass.mesh import QuasiUniformSphericalMeshStep -# from compass.mesh import QuasiUniformSphericalMeshStep - -class WavesUostFiles(): +class WavesUostFiles(QuasiUniformSphericalMeshStep): """ A step for creating the unresolved obstacles file for wave mesh """ - # pass + def __init__(self, test_case, ocean_mesh, name='uost_files', subdir=None): - def __init__(self, test_case): - """ - Create a new step - """ - super().__init__(test_case, name='uost_files') + super().__init__(test_case=test_case, name=name, subdir=subdir, + cell_width=None) + + # other things INIT should do? + base_mesh_path = ocean_mesh.steps['base_mesh'].path + self.add_input_file( + filename='ocean_mesh.nc', + work_dir_target=f'{base_mesh_path}/base_mesh.nc') + + self.add_input_file( + filename='etopo1_180.nc', + target='etopo1_180.nc', + database='bathymetry_database') def run(self): """ - Run this step + Create unresolved obstacles for wave mesh and spectral resolution """ - super().run() dirs = np.linspace(0, 2 * np.pi, 36) nfreq = 50 # ET NOTE: this should be flexible @@ -37,7 +43,9 @@ def run(self): frqfactor = 1.10 elif (nfreq == 25): frqfactor = 1.147 - # PRINT SOME ERROR MESSAGE FOR NON-supoorted spectral reslution + else: + print("ERROR: Spectral resolution not supported.") + print("Number of wave freqencies must be 25, 36, or 50.") freqs = [minfrq * (frqfactor ** i) for i in range(1, nfreq + 1)] @@ -69,32 +77,3 @@ def run(self): abEstimateAndSaveTriangularEtopo1( dirs, freqs, gridname, triMeshSpec, etopoFilePath, outputDestDir, nParWorker, abOptions=opt) - - # def build_cell_width_lat_lon(self): - # """ - # Create cell width array for this mesh on a regular latitude-longitude - # grid - - # Returns - # ------- - # cellWidth : numpy.array - # m x n array of cell width in km - - # lon : numpy.array - # longitude in degrees (length n and between -180 and 180) - - # lat : numpy.array - # longitude in degrees (length m and between -90 and 90) - # """ - - # dlon = 10. - # dlat = 0.1 - # nlon = int(360. / dlon) + 1 - # nlat = int(180. / dlat) + 1 - # lon = np.linspace(-180., 180., nlon) - # lat = np.linspace(-90., 90., nlat) - - # cellWidthVsLat = mdt.EC_CellWidthVsLat(lat) - # cellWidth = np.outer(cellWidthVsLat, np.ones([1, lon.size])) - - # return cellWidth, lon, lat From 012c86b8731079e167b4c48c68cc4cf1304cf8f1 Mon Sep 17 00:00:00 2001 From: Erin Thomas Date: Wed, 17 Apr 2024 15:45:09 -0500 Subject: [PATCH 06/26] fix to __init__ in call to uost --- compass/ocean/tests/global_ocean/wave_mesh/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 0d550c8cc4..d04878019f 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -37,8 +37,9 @@ def __init__(self, test_group, ocean_mesh): base_mesh_step = WavesBaseMesh(test_case=self, ocean_mesh=ocean_mesh) self.add_step(base_mesh_step) - - uost_file_step = WavesUostFiles(test_case=self) + + uost_file_step = WavesUostFiles(test_case=self, + ocean_mesh=ocean_mesh) self.add_step(uost_file_step) def configure(self, config=None): From cc54225c559ae08be7a9067cd50e391722366d3c Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Wed, 17 Apr 2024 15:31:20 -0500 Subject: [PATCH 07/26] Add cull mesh step --- .../tests/global_ocean/wave_mesh/__init__.py | 6 +++ .../tests/global_ocean/wave_mesh/cull_mesh.py | 42 +++++++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index d04878019f..2f788fa400 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -5,6 +5,7 @@ from compass.ocean.tests.global_ocean.wave_mesh.uost_files import ( WavesUostFiles, ) +from compass.ocean.tests.global_ocean.wave_mesh.cull_mesh import WavesCullMesh from compass.testcase import TestCase # from compass.validate import compare_variables @@ -38,6 +39,11 @@ def __init__(self, test_group, ocean_mesh): ocean_mesh=ocean_mesh) self.add_step(base_mesh_step) + cull_mesh_step = WavesCullMesh(test_case=self, + ocean_mesh=ocean_mesh, + wave_base_mesh=base_mesh_step) + self.add_step(cull_mesh_step) + uost_file_step = WavesUostFiles(test_case=self, ocean_mesh=ocean_mesh) self.add_step(uost_file_step) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py new file mode 100644 index 0000000000..c6f01e820e --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py @@ -0,0 +1,42 @@ +from mpas_tools.logging import check_call + +from compass import Step + + +class WavesCullMesh(Step): + """ + A step for creating wave mesh based on an ocean mesh + """ + def __init__(self, test_case, ocean_mesh, wave_base_mesh, + name='cull_mesh', subdir=None): + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + culled_mesh_path = ocean_mesh.steps['cull_mesh'].path + self.add_input_file( + filename='ocean_culled_mesh.nc', + work_dir_target=f'{culled_mesh_path}/culled_mesh.nc') + + wave_base_mesh_path = wave_base_mesh.path + self.add_input_file( + filename='wave_base_mesh.nc', + work_dir_target=f'{wave_base_mesh_path}/wave_base_mesh.nc') + + def setup(self): + + super().setup() + + f = open(f'{self.work_dir}/cull_waves_mesh.nml', 'w') + f.write('&inputs\n') + f.write(" waves_mesh_file = 'wave_base_mesh.nc'\n") + f.write(" ocean_mesh_file = 'ocean_culled_mesh.nc'\n") + f.write("/\n") + f.write("&output\n") + f.write(" waves_mesh_culled_vtk = 'wave_mesh_culled.vtk'\n") + f.write(" waves_mesh_culled_gmsh = 'wave_mesh_culled.msh'\n") + f.write("/\n") + f.close() + + def run(self): + + check_call('ocean_cull_wave_mesh', logger=self.logger) From 115f539de41598dfd41f4bc9e377a990d46b2973 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Wed, 17 Apr 2024 16:20:31 -0500 Subject: [PATCH 08/26] Add rotate step --- .../tests/global_ocean/wave_mesh/__init__.py | 9 ++++- .../global_ocean/wave_mesh/rotate_mesh.py | 36 +++++++++++++++++++ 2 files changed, 44 insertions(+), 1 deletion(-) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 2f788fa400..f5e89ad243 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -2,10 +2,13 @@ # get_author_and_email_from_git, # ) from compass.ocean.tests.global_ocean.wave_mesh.base_mesh import WavesBaseMesh +from compass.ocean.tests.global_ocean.wave_mesh.cull_mesh import WavesCullMesh +from compass.ocean.tests.global_ocean.wave_mesh.rotate_mesh import ( + WavesRotateMesh, +) from compass.ocean.tests.global_ocean.wave_mesh.uost_files import ( WavesUostFiles, ) -from compass.ocean.tests.global_ocean.wave_mesh.cull_mesh import WavesCullMesh from compass.testcase import TestCase # from compass.validate import compare_variables @@ -44,6 +47,10 @@ def __init__(self, test_group, ocean_mesh): wave_base_mesh=base_mesh_step) self.add_step(cull_mesh_step) + rotate_mesh_step = WavesRotateMesh(test_case=self, + wave_culled_mesh=cull_mesh_step) + self.add_step(rotate_mesh_step) + uost_file_step = WavesUostFiles(test_case=self, ocean_mesh=ocean_mesh) self.add_step(uost_file_step) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py new file mode 100644 index 0000000000..cb1862204f --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py @@ -0,0 +1,36 @@ +from mpas_tools.logging import check_call + +from compass import Step + + +class WavesRotateMesh(Step): + """ + A step for creating wave mesh based on an ocean mesh + """ + def __init__(self, test_case, wave_culled_mesh, + name='rotate_mesh', subdir=None): + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + wave_culled_mesh_path = wave_culled_mesh.path + self.add_input_file( + filename='wave_mesh_culled.msh', + work_dir_target=f'{wave_culled_mesh_path}/wave_mesh_culled.msh') + + def setup(self): + + super().setup() + + f = open(f'{self.work_dir}/rotate.nml', 'w') + f.write('&inputs\n') + f.write(" LON_POLE = -42.8906d0\n") + f.write(" LAT_POLE = 72.3200d0\n") + f.write(" wind_file = 'null'\n") + f.write(" mesh_file = 'wave_mesh_culled.msh'\n") + f.write(" mesh_file_out = 'wave_culled_mesh_RTD.msh'\n") + f.write("/") + f.close() + + def run(self): + + check_call('ocean_rotate_wave_mesh', logger=self.logger) From b062dd58723b7ec2dc585168c507890eddd8b4c2 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Wed, 17 Apr 2024 16:40:41 -0500 Subject: [PATCH 09/26] Minor changes to uost step --- .../tests/global_ocean/wave_mesh/__init__.py | 2 +- .../global_ocean/wave_mesh/uost_files.py | 24 +++++++++---------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index f5e89ad243..09adc6d939 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -52,7 +52,7 @@ def __init__(self, test_group, ocean_mesh): self.add_step(rotate_mesh_step) uost_file_step = WavesUostFiles(test_case=self, - ocean_mesh=ocean_mesh) + wave_culled_mesh=cull_mesh_step) self.add_step(uost_file_step) def configure(self, config=None): diff --git a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py index e0e99bbacc..1e49cead51 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py @@ -6,23 +6,23 @@ ) # importing from alphaBetaLab the needed components from alphaBetaLab.alphaBetaLab.abOptionManager import abOptions -from compass.mesh import QuasiUniformSphericalMeshStep +from compass import Step -class WavesUostFiles(QuasiUniformSphericalMeshStep): +class WavesUostFiles(Step): """ A step for creating the unresolved obstacles file for wave mesh """ - def __init__(self, test_case, ocean_mesh, name='uost_files', subdir=None): + def __init__(self, test_case, wave_culled_mesh, + name='uost_files', subdir=None): - super().__init__(test_case=test_case, name=name, subdir=subdir, - cell_width=None) + super().__init__(test_case=test_case, name=name, subdir=subdir) # other things INIT should do? - base_mesh_path = ocean_mesh.steps['base_mesh'].path + culled_mesh_path = wave_culled_mesh.path self.add_input_file( - filename='ocean_mesh.nc', - work_dir_target=f'{base_mesh_path}/base_mesh.nc') + filename='wave_mesh_culled.msh', + work_dir_target=f'{culled_mesh_path}/wave_mesh_culled.msh') self.add_input_file( filename='etopo1_180.nc', @@ -50,17 +50,17 @@ def run(self): freqs = [minfrq * (frqfactor ** i) for i in range(1, nfreq + 1)] # definition of the spatial mesh - gridname = 'glo_unst' - mshfile = 'global.msh' + gridname = 'glo_unst' # 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' + etopoFilePath = 'etopo1_180.nc' # output directory - outputDestDir = './output/' + outputDestDir = 'output/' # number of cores for parallel computing nParWorker = 1 From dadb9fb9614484f1e8e99dc89071943d7ddde71a Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Wed, 17 Apr 2024 20:18:13 -0500 Subject: [PATCH 10/26] Use initial_state.nc as ocean mesh --- compass/ocean/tests/global_ocean/__init__.py | 3 ++- compass/ocean/tests/global_ocean/wave_mesh/__init__.py | 6 +++--- compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py | 8 ++++---- compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py | 6 +++--- compass/ocean/tests/global_ocean/wave_mesh/uost_files.py | 2 +- 5 files changed, 13 insertions(+), 12 deletions(-) diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index f11572bc89..6393d7cd4f 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -175,4 +175,5 @@ def _add_tests(self, mesh_names, high_res_topography=True, if include_wave_mesh: self.add_test_case( - WaveMesh(test_group=self, ocean_mesh=mesh_test)) + WaveMesh(test_group=self, ocean_mesh=mesh_test, + ocean_init=init_test)) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 09adc6d939..61bc079c30 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -21,7 +21,7 @@ class WaveMesh(TestCase): Attributes ---------- """ - def __init__(self, test_group, ocean_mesh): + def __init__(self, test_group, ocean_mesh, ocean_init): """ Create test case for creating a global MPAS-Ocean mesh @@ -39,11 +39,11 @@ def __init__(self, test_group, ocean_mesh): self.mesh_config_filename = 'wave_mesh.cfg' base_mesh_step = WavesBaseMesh(test_case=self, - ocean_mesh=ocean_mesh) + ocean_mesh=ocean_init) self.add_step(base_mesh_step) cull_mesh_step = WavesCullMesh(test_case=self, - ocean_mesh=ocean_mesh, + ocean_mesh=ocean_init, wave_base_mesh=base_mesh_step) self.add_step(cull_mesh_step) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index 2bdc994b04..c5a00c83f4 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -23,10 +23,10 @@ def __init__(self, test_case, ocean_mesh, name='base_mesh', subdir=None): super().__init__(test_case=test_case, name=name, subdir=subdir, cell_width=None) - base_mesh_path = ocean_mesh.steps['base_mesh'].path + mesh_path = ocean_mesh.steps['initial_state'].path self.add_input_file( filename='ocean_mesh.nc', - work_dir_target=f'{base_mesh_path}/base_mesh.nc') + work_dir_target=f'{mesh_path}/initial_state.nc') def setup(self): @@ -182,8 +182,8 @@ def limit_spacing_gradient(self, lon, lat, cell_width, dhdx): spac = jigsawpy.jigsaw_msh_t() spac.mshID = "ellipsoid-grid" - spac.xgrid = lon - spac.ygrid = lat + spac.xgrid = np.radians(lon) + spac.ygrid = np.radians(lat) spac.value = cell_width spac.slope = np.full(spac.value.shape, dhdx) jigsawpy.savemsh(opts.hfun_file, spac) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py index c6f01e820e..1994b9054b 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py @@ -12,15 +12,15 @@ def __init__(self, test_case, ocean_mesh, wave_base_mesh, super().__init__(test_case=test_case, name=name, subdir=subdir) - culled_mesh_path = ocean_mesh.steps['cull_mesh'].path + culled_mesh_path = ocean_mesh.steps['initial_state'].path self.add_input_file( filename='ocean_culled_mesh.nc', - work_dir_target=f'{culled_mesh_path}/culled_mesh.nc') + work_dir_target=f'{culled_mesh_path}/initial_state.nc') wave_base_mesh_path = wave_base_mesh.path self.add_input_file( filename='wave_base_mesh.nc', - work_dir_target=f'{wave_base_mesh_path}/wave_base_mesh.nc') + work_dir_target=f'{wave_base_mesh_path}/base_mesh.nc') def setup(self): diff --git a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py index 1e49cead51..76ddb5c357 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py @@ -63,7 +63,7 @@ def run(self): outputDestDir = 'output/' # number of cores for parallel computing - nParWorker = 1 + nParWorker = 1 # SB NOTE: we should set this to the of cores on a node # this option indicates that the computation # should be skipped for cells smaller than 3 km From 6f6c14e2b0e0dbedad378955791900a94e2b9d15 Mon Sep 17 00:00:00 2001 From: Erin Thomas Date: Wed, 17 Apr 2024 21:58:47 -0500 Subject: [PATCH 11/26] scrip file creation --- .../global_ocean/wave_mesh/scrip_file.py | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py new file mode 100644 index 0000000000..31b7476c9c --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -0,0 +1,37 @@ +from mpas_tools.logging import check_call + +from compass import Step + + +class WavesScripFile(Step): + """ + A step for creating the scrip file for the wave mesh + """ + def __init__(self, test_case, wave_culled_mesh, + name='scrip_file', subdir=None): + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + wave_culled_mesh_path = wave_culled_mesh.path + self.add_input_file( + filename='wave_mesh_culled.msh', + work_dir_target=f'{wave_culled_mesh_path}/wave_mesh_culled.msh') + + def setup(self): + + super().setup() + + f = open(f'{self.work_dir}/scrip.nml', 'w') + f.write('&inputs\n') + f.write(" waves_mesh_file = 'waves_mesh_culled.msh'\n") + f.write("/") + f.write('&outputs\n') + f.write(" waves_scrip_file = 'waves_mesh_scrip.nc'\n") + f.write("/") + f.close() + + def run(self): + """ + Create scrip files for wave mesh + """ + check_call('ocean_scrip_wave_mesh', logger=self.logger) From 30bc773e7148ede94c2d75fffe7a786feda29a45 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Thu, 18 Apr 2024 11:13:47 -0500 Subject: [PATCH 12/26] Fix for scrip step --- compass/ocean/tests/global_ocean/wave_mesh/__init__.py | 7 +++++++ compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py | 6 +++--- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 61bc079c30..15a9811350 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -6,6 +6,9 @@ from compass.ocean.tests.global_ocean.wave_mesh.rotate_mesh import ( WavesRotateMesh, ) +from compass.ocean.tests.global_ocean.wave_mesh.scrip_file import ( + WavesScripFile, +) from compass.ocean.tests.global_ocean.wave_mesh.uost_files import ( WavesUostFiles, ) @@ -47,6 +50,10 @@ def __init__(self, test_group, ocean_mesh, ocean_init): wave_base_mesh=base_mesh_step) self.add_step(cull_mesh_step) + scrip_file_step = WavesScripFile(test_case=self, + wave_culled_mesh=cull_mesh_step) + self.add_step(scrip_file_step) + rotate_mesh_step = WavesRotateMesh(test_case=self, wave_culled_mesh=cull_mesh_step) self.add_step(rotate_mesh_step) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index 31b7476c9c..468a15af2b 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -23,10 +23,10 @@ def setup(self): f = open(f'{self.work_dir}/scrip.nml', 'w') f.write('&inputs\n') - f.write(" waves_mesh_file = 'waves_mesh_culled.msh'\n") - f.write("/") + f.write(" waves_mesh_file = 'wave_mesh_culled.msh'\n") + f.write("/\n") f.write('&outputs\n') - f.write(" waves_scrip_file = 'waves_mesh_scrip.nc'\n") + f.write(" waves_scrip_file = 'wave_mesh_scrip.nc'\n") f.write("/") f.close() From b72581e207ec690e6e0f086f564a8f79b217dd9b Mon Sep 17 00:00:00 2001 From: Erin Thomas Date: Thu, 18 Apr 2024 13:44:22 -0500 Subject: [PATCH 13/26] initial remap step --- .../tests/global_ocean/wave_mesh/__init__.py | 9 ++ .../global_ocean/wave_mesh/remap_files.py | 106 ++++++++++++++++++ 2 files changed, 115 insertions(+) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/remap_files.py diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 15a9811350..ac055a447e 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -1,8 +1,12 @@ # from compass.ocean.tests.global_ocean.metadata import ( # get_author_and_email_from_git, # ) +from compass.ocean.tests.global_ocean.files_for_e3sm.scrip import Scrip from compass.ocean.tests.global_ocean.wave_mesh.base_mesh import WavesBaseMesh from compass.ocean.tests.global_ocean.wave_mesh.cull_mesh import WavesCullMesh +from compass.ocean.tests.global_ocean.wave_mesh.remap_files import ( + WavesRemapFiles, +) from compass.ocean.tests.global_ocean.wave_mesh.rotate_mesh import ( WavesRotateMesh, ) @@ -54,6 +58,11 @@ def __init__(self, test_group, ocean_mesh, ocean_init): wave_culled_mesh=cull_mesh_step) self.add_step(scrip_file_step) + remap_file_step = WavesRemapFiles(test_case=self, + wave_scrip=scrip_file_step, + ocean_scrip=Scrip) + self.add_step(remap_file_step) + rotate_mesh_step = WavesRotateMesh(test_case=self, wave_culled_mesh=cull_mesh_step) self.add_step(rotate_mesh_step) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py b/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py new file mode 100644 index 0000000000..ff391a00d0 --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py @@ -0,0 +1,106 @@ +import os +import pprint +import subprocess + +import yaml + +from compass import Step + + +class WavesRemapFiles(Step): + """ + A step for creating remapping files for wave mesh + """ + def __init__(self, test_case, wave_scrip, ocean_scrip, + name='remap_files', subdir=None): + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + wave_scrip_file_path = wave_scrip.path + self.add_input_file( + filename='wave_scrip.nc', + work_dir_target=f'{wave_scrip_file_path}/wave_scrip.nc') + + ocean_scrip_file_path = ocean_scrip.path + self.add_input_file( + filename='ocean_scrip.nc', + work_dir_target=f'{ocean_scrip_file_path}/ocean_scrip.nc') + + def setup(self): + + super().setup() + # TO DO: mesh names should be flexible not hard-coded. + f = open(f'{self.work_dir}/make_remapping_files.config', 'w') + f.write("grid1 : 'scrip.nc'\n") + f.write("grid1_shortname : 'wQU225IcoswISC30E3r5'\n") + f.write("reg1 : True\n") + f.write("\n") + f.write("grid2 : 'ocean.IcoswISC30E3r5.mask.scrip.20231120.nc'\n") + f.write("grid2_shortname : 'IcoswISC30E3r5'\n") + f.write("reg2 : True\n") + f.write("\n") + f.write("datestamp : 20240411\n") + f.write("\n") + f.write("map_types : ['conserve','bilinear','neareststod']\n") + f.close() + + def make_remapping_files(grid1, grid2, + grid1_shortname, grid2_shortname, + datestamp, reg1, reg2, map_type, nprocs=1): + if map_type == 'conserve': + map_abbrev = 'aave' + elif map_type == 'bilinear': + map_abbrev = 'blin' + elif map_type == 'patch': + map_abbrev = 'patc' + elif map_type == 'neareststod': + map_abbrev = 'nstod' + elif map_type == 'nearsetdtos': + map_abbrev = 'ndtos' + else: + print('map type not recognized') + raise SystemExit(0) + + flags = ' --ignore_unmapped --ignore_degenerate' + if reg1: + flags += ' --src_regional' + if reg2: + flags += ' --dst_regional' + + map_name = 'map_' + grid1_shortname + '_TO_' + grid2_shortname +\ + '_' + map_abbrev + '.' + str(datestamp) + '.nc' + subprocess.call('srun -n ' + str(nprocs) + + ' ESMF_RegridWeightGen --source ' + grid1 + + ' --destination ' + grid2 + + ' --method ' + map_type + + ' --weight ' + map_name + + flags, shell=True) + + flags = ' --ignore_unmapped --ignore_degenerate' + if reg1: + flags += ' --dst_regional' + if reg2: + flags += ' --src_regional' + + map_name = 'map_' + grid2_shortname + '_TO_' + grid1_shortname +\ + '_' + map_abbrev + '.' + str(datestamp) + '.nc' + subprocess.call('srun -n ' + str(nprocs) + + ' ESMF_RegridWeightGen --source ' + grid2 + + ' --destination ' + grid1 + + ' --method ' + map_type + + ' --weight ' + map_name + + flags, shell=True) + + def run(self): + pwd = os.getcwd() + f = open(pwd + '/make_remapping_files.config') + cfg = yaml.load(f, Loader=yaml.Loader) + pprint.pprint(cfg) + for map_type in cfg['map_types']: + self.make_remapping_files(cfg['grid1'], + cfg['grid2'], + cfg['grid1_shortname'], + cfg['grid2_shortname'], + cfg['datestamp'], + cfg['reg1'], cfg['reg2'], + map_type) From 08f9d12c61f86036a28d3abee6c0f87e49f9e5c0 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Thu, 18 Apr 2024 15:59:25 -0500 Subject: [PATCH 14/26] Namelist updates --- .../ocean/tests/global_ocean/wave_mesh/rotate_mesh.py | 10 +++++----- .../ocean/tests/global_ocean/wave_mesh/scrip_file.py | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py index cb1862204f..0388f40d7a 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py @@ -23,11 +23,11 @@ def setup(self): f = open(f'{self.work_dir}/rotate.nml', 'w') f.write('&inputs\n') - f.write(" LON_POLE = -42.8906d0\n") - f.write(" LAT_POLE = 72.3200d0\n") - f.write(" wind_file = 'null'\n") - f.write(" mesh_file = 'wave_mesh_culled.msh'\n") - f.write(" mesh_file_out = 'wave_culled_mesh_RTD.msh'\n") + f.write(" LON_POLE = -42.8906d0\n") + f.write(" LAT_POLE = 72.3200d0\n") + f.write(" wind_file = 'null'\n") + f.write(" mesh_file = 'wave_mesh_culled.msh'\n") + f.write(" mesh_file_out = 'wave_culled_mesh_RTD.msh'\n") f.write("/") f.close() diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index 468a15af2b..332751a23e 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -23,10 +23,10 @@ def setup(self): f = open(f'{self.work_dir}/scrip.nml', 'w') f.write('&inputs\n') - f.write(" waves_mesh_file = 'wave_mesh_culled.msh'\n") + f.write(" waves_mesh_file = 'wave_mesh_culled.msh'\n") f.write("/\n") f.write('&outputs\n') - f.write(" waves_scrip_file = 'wave_mesh_scrip.nc'\n") + f.write(" waves_scrip_file = 'wave_mesh_scrip.nc'\n") f.write("/") f.close() From af105551b38b7eb941951fb3a471c8ae753af631 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 19 Apr 2024 11:41:29 -0500 Subject: [PATCH 15/26] Add new lines at end of rotate and scrip nml files - Avoids end of file error in namelist read --- compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py | 2 +- compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py index 0388f40d7a..e72d254227 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py @@ -28,7 +28,7 @@ def setup(self): f.write(" wind_file = 'null'\n") f.write(" mesh_file = 'wave_mesh_culled.msh'\n") f.write(" mesh_file_out = 'wave_culled_mesh_RTD.msh'\n") - f.write("/") + f.write("/\n") f.close() def run(self): diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index 332751a23e..0e60ccd98f 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -27,7 +27,7 @@ def setup(self): f.write("/\n") f.write('&outputs\n') f.write(" waves_scrip_file = 'wave_mesh_scrip.nc'\n") - f.write("/") + f.write("/\n") f.close() def run(self): From b86b6f435ab01708ae7223a8e593a1071f3b4ae9 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 19 Apr 2024 14:21:28 -0500 Subject: [PATCH 16/26] Update remap step --- compass/ocean/tests/global_ocean/__init__.py | 10 +-- .../tests/global_ocean/wave_mesh/__init__.py | 5 +- .../global_ocean/wave_mesh/remap_files.py | 84 +++++++------------ 3 files changed, 38 insertions(+), 61 deletions(-) diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index 6393d7cd4f..adbb797f41 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -121,10 +121,10 @@ def _add_tests(self, mesh_names, high_res_topography=True, time_integrator=time_integrator) self.add_test_case(dynamic_adjustment_test) - self.add_test_case( - FilesForE3SM( - test_group=self, mesh=mesh_test, init=init_test, - dynamic_adjustment=dynamic_adjustment_test)) + e3sm_files_test = FilesForE3SM( + test_group=self, mesh=mesh_test, init=init_test, + dynamic_adjustment=dynamic_adjustment_test) + self.add_test_case(e3sm_files_test) if include_rk4: time_integrator = 'RK4' @@ -176,4 +176,4 @@ def _add_tests(self, mesh_names, high_res_topography=True, if include_wave_mesh: self.add_test_case( WaveMesh(test_group=self, ocean_mesh=mesh_test, - ocean_init=init_test)) + ocean_init=init_test, ocean_e3sm=e3sm_files_test)) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index ac055a447e..d77b7a10e0 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -1,7 +1,6 @@ # from compass.ocean.tests.global_ocean.metadata import ( # get_author_and_email_from_git, # ) -from compass.ocean.tests.global_ocean.files_for_e3sm.scrip import Scrip from compass.ocean.tests.global_ocean.wave_mesh.base_mesh import WavesBaseMesh from compass.ocean.tests.global_ocean.wave_mesh.cull_mesh import WavesCullMesh from compass.ocean.tests.global_ocean.wave_mesh.remap_files import ( @@ -28,7 +27,7 @@ class WaveMesh(TestCase): Attributes ---------- """ - def __init__(self, test_group, ocean_mesh, ocean_init): + def __init__(self, test_group, ocean_mesh, ocean_init, ocean_e3sm): """ Create test case for creating a global MPAS-Ocean mesh @@ -60,7 +59,7 @@ def __init__(self, test_group, ocean_mesh, ocean_init): remap_file_step = WavesRemapFiles(test_case=self, wave_scrip=scrip_file_step, - ocean_scrip=Scrip) + ocean_e3sm=ocean_e3sm) self.add_step(remap_file_step) rotate_mesh_step = WavesRotateMesh(test_case=self, diff --git a/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py b/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py index ff391a00d0..64afe1fbac 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py @@ -1,8 +1,6 @@ -import os -import pprint -import subprocess +from datetime import date -import yaml +from mpas_tools.logging import check_call from compass import Step @@ -11,7 +9,7 @@ class WavesRemapFiles(Step): """ A step for creating remapping files for wave mesh """ - def __init__(self, test_case, wave_scrip, ocean_scrip, + def __init__(self, test_case, wave_scrip, ocean_e3sm, name='remap_files', subdir=None): super().__init__(test_case=test_case, name=name, subdir=subdir) @@ -19,32 +17,14 @@ def __init__(self, test_case, wave_scrip, ocean_scrip, wave_scrip_file_path = wave_scrip.path self.add_input_file( filename='wave_scrip.nc', - work_dir_target=f'{wave_scrip_file_path}/wave_scrip.nc') + work_dir_target=f'{wave_scrip_file_path}/wave_mesh_scrip.nc') - ocean_scrip_file_path = ocean_scrip.path + ocean_scrip_file_path = ocean_e3sm.steps['scrip'].path self.add_input_file( filename='ocean_scrip.nc', work_dir_target=f'{ocean_scrip_file_path}/ocean_scrip.nc') - def setup(self): - - super().setup() - # TO DO: mesh names should be flexible not hard-coded. - f = open(f'{self.work_dir}/make_remapping_files.config', 'w') - f.write("grid1 : 'scrip.nc'\n") - f.write("grid1_shortname : 'wQU225IcoswISC30E3r5'\n") - f.write("reg1 : True\n") - f.write("\n") - f.write("grid2 : 'ocean.IcoswISC30E3r5.mask.scrip.20231120.nc'\n") - f.write("grid2_shortname : 'IcoswISC30E3r5'\n") - f.write("reg2 : True\n") - f.write("\n") - f.write("datestamp : 20240411\n") - f.write("\n") - f.write("map_types : ['conserve','bilinear','neareststod']\n") - f.close() - - def make_remapping_files(grid1, grid2, + def make_remapping_files(self, grid1, grid2, grid1_shortname, grid2_shortname, datestamp, reg1, reg2, map_type, nprocs=1): if map_type == 'conserve': @@ -67,14 +47,14 @@ def make_remapping_files(grid1, grid2, if reg2: flags += ' --dst_regional' - map_name = 'map_' + grid1_shortname + '_TO_' + grid2_shortname +\ - '_' + map_abbrev + '.' + str(datestamp) + '.nc' - subprocess.call('srun -n ' + str(nprocs) + - ' ESMF_RegridWeightGen --source ' + grid1 + - ' --destination ' + grid2 + - ' --method ' + map_type + - ' --weight ' + map_name + - flags, shell=True) + map_name = f'map_{grid1_shortname}_TO_{grid2_shortname}'\ + f'_{map_abbrev}.{datestamp}.nc' + check_call(f'srun -n {nprocs}' + f' ESMF_RegridWeightGen --source {grid1}' + f' --destination {grid2}' + f' --method {map_type}' + f' --weight {map_name} {flags}', + logger=self.logger) flags = ' --ignore_unmapped --ignore_degenerate' if reg1: @@ -82,25 +62,23 @@ def make_remapping_files(grid1, grid2, if reg2: flags += ' --src_regional' - map_name = 'map_' + grid2_shortname + '_TO_' + grid1_shortname +\ - '_' + map_abbrev + '.' + str(datestamp) + '.nc' - subprocess.call('srun -n ' + str(nprocs) + - ' ESMF_RegridWeightGen --source ' + grid2 + - ' --destination ' + grid1 + - ' --method ' + map_type + - ' --weight ' + map_name + - flags, shell=True) + map_name = f'map_{grid2_shortname}_TO_{grid1_shortname}'\ + f'_{map_abbrev}.{datestamp}.nc' + check_call(f'srun -n {nprocs}' + f' ESMF_RegridWeightGen --source {grid2}' + f' --destination {grid1}' + f' --method {map_type}' + f' --weight {map_name} {flags}', + logger=self.logger) def run(self): - pwd = os.getcwd() - f = open(pwd + '/make_remapping_files.config') - cfg = yaml.load(f, Loader=yaml.Loader) - pprint.pprint(cfg) - for map_type in cfg['map_types']: - self.make_remapping_files(cfg['grid1'], - cfg['grid2'], - cfg['grid1_shortname'], - cfg['grid2_shortname'], - cfg['datestamp'], - cfg['reg1'], cfg['reg2'], + today = date.today() + creation_date = today.strftime("%Y%m%d") + for map_type in ['conserve', 'bilinear', 'neareststod']: + self.make_remapping_files('wave_scrip.nc', + 'ocean_scrip.nc', + 'wave', # This should be more specific + 'ocean', # This should be more spedific + creation_date, + True, True, map_type) From 504d69adb792da746bdac037922d0ebcb44ccae2 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 19 Apr 2024 14:22:37 -0500 Subject: [PATCH 17/26] Fix resolution specification in base_mesh step --- .../tests/global_ocean/wave_mesh/base_mesh.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index c5a00c83f4..fe09a8b5fd 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -63,6 +63,7 @@ def build_cell_width_lat_lon(self): ylat = np.linspace(lat_min, lat_max, nlat) earth_radius = constants['SHR_CONST_REARTH'] / km + # SB NOTE: need to use GHSSH from cartopy shapefiles = ['GSHHS_l_L1.shp', 'GSHHS_l_L6.shp'] cell_width = self.cell_widthVsLatLon(xlon, ylat, shapefiles, earth_radius, 'ocean_mesh.nc') @@ -71,6 +72,10 @@ def build_cell_width_lat_lon(self): self.create_initial_points('ocean_mesh.nc', xlon, ylat, cell_width, earth_radius, self.opts.init_file) + hfun_slope_lim = self.config.getfloat('wave_mesh', 'hfun_slope_lim') + cell_width = self.limit_spacing_gradient(xlon, ylat, cell_width, + earth_radius, hfun_slope_lim) + return cell_width, xlon, ylat def cell_widthVsLatLon(self, lon, lat, coast_shapefiles, @@ -140,9 +145,6 @@ def cell_widthVsLatLon(self, lon, lat, coast_shapefiles, (D < dist_threshold_global)) cell_width[idxx, idxy] = hfun_grd[idxx, idxy] - hfun_slope_lim = config.getfloat('wave_mesh', 'hfun_slope_lim') - self.limit_spacing_gradient(lon, lat, cell_width, hfun_slope_lim) - # Plot fig = plt.figure(figsize=(16, 8)) @@ -171,7 +173,8 @@ def cell_widthVsLatLon(self, lon, lat, coast_shapefiles, return cell_width - def limit_spacing_gradient(self, lon, lat, cell_width, dhdx): + def limit_spacing_gradient(self, lon, lat, cell_width, + sphere_radius, dhdx): print("Smoothing h(x) via |dh/dx| limits...") @@ -182,10 +185,11 @@ def limit_spacing_gradient(self, lon, lat, cell_width, dhdx): spac = jigsawpy.jigsaw_msh_t() spac.mshID = "ellipsoid-grid" + spac.radii = np.full(3, sphere_radius, dtype=spac.REALS_t) spac.xgrid = np.radians(lon) spac.ygrid = np.radians(lat) - spac.value = cell_width - spac.slope = np.full(spac.value.shape, dhdx) + spac.value = cell_width.astype(spac.REALS_t) + spac.slope = np.full(spac.value.shape, dhdx, dtype=spac.REALS_t) jigsawpy.savemsh(opts.hfun_file, spac) jigsawpy.cmd.marche(opts, spac) @@ -306,7 +310,6 @@ def create_initial_points(self, meshfile, lon, lat, hfunction, lat = np.append(lat, 0.5 * np.pi) npt = lon.size - print(npt) # Change to Cartesian coordinates x, y, z = self.lonlat2xyz(lon, lat, sphere_radius) From 68476150e5b1b0be959a9bc3e6ffc04eed37df25 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 19 Apr 2024 16:04:44 -0500 Subject: [PATCH 18/26] Use GSHHS coastlines through cartopy --- .../tests/global_ocean/wave_mesh/base_mesh.py | 38 +++++++------------ 1 file changed, 13 insertions(+), 25 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index fe09a8b5fd..726df345b0 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -1,12 +1,10 @@ import timeit -import cartopy # noqa: F401 -import cartopy.crs as ccrs # noqa: F401 -import cartopy.feature as cfeature # noqa: F401 +import cartopy.crs as ccrs +import cartopy.feature as cfeature import jigsawpy import matplotlib.pyplot as plt import numpy as np -import shapefile from mpas_tools.cime.constants import constants from netCDF4 import Dataset from scipy import interpolate, spatial @@ -63,9 +61,7 @@ def build_cell_width_lat_lon(self): ylat = np.linspace(lat_min, lat_max, nlat) earth_radius = constants['SHR_CONST_REARTH'] / km - # SB NOTE: need to use GHSSH from cartopy - shapefiles = ['GSHHS_l_L1.shp', 'GSHHS_l_L6.shp'] - cell_width = self.cell_widthVsLatLon(xlon, ylat, shapefiles, + cell_width = self.cell_widthVsLatLon(xlon, ylat, earth_radius, 'ocean_mesh.nc') cell_width = cell_width / km @@ -78,8 +74,7 @@ def build_cell_width_lat_lon(self): return cell_width, xlon, ylat - def cell_widthVsLatLon(self, lon, lat, coast_shapefiles, - sphere_radius, ocean_mesh): + def cell_widthVsLatLon(self, lon, lat, sphere_radius, ocean_mesh): config = self.config @@ -129,7 +124,7 @@ def cell_widthVsLatLon(self, lon, lat, coast_shapefiles, bathy_grd = np.reshape(bathy_interp, (nlat, nlon)) # Get distance to coasts - D = self.distance_to_shapefile_points(coast_shapefiles, lon, lat, + D = self.distance_to_shapefile_points(lon, lat, sphere_radius, reggrid=True) # Apply refined region criteria @@ -196,27 +191,20 @@ def limit_spacing_gradient(self, lon, lat, cell_width, return spac.value - def distance_to_shapefile_points(self, shpfiles, lon, lat, + def distance_to_shapefile_points(self, lon, lat, sphere_radius, reggrid=False): # Get coastline coordinates from shapefile - pt_list = [] - for shpfile in shpfiles: - - sf = shapefile.Reader(shpfile) - shapes = sf.shapes() - # records = sf.records() - n = len(sf.shapes()) + features = cfeature.GSHHSFeature(scale='l', levels=[1, 6]) - for i in range(n): - - points = shapes[i].points + pt_list = [] + for feature in features.geometries(): - if len(points) < 20: - continue + if feature.length < 20: + continue - for pt in points: - pt_list.append([pt[0], pt[1]]) + for coord in feature.exterior.coords: + pt_list.append([coord[0], coord[1]]) coast_pts = np.radians(np.array(pt_list)) From 8a8c9f56f7a0eea60bd6a0b3db988699206d066c Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 19 Apr 2024 16:45:12 -0500 Subject: [PATCH 19/26] Use jinja templates for writing namelists --- .../tests/global_ocean/wave_mesh/cull_mesh.py | 22 ++++++++++--------- .../global_ocean/wave_mesh/cull_mesh.template | 8 +++++++ .../global_ocean/wave_mesh/rotate_mesh.py | 21 ++++++++++-------- .../wave_mesh/rotate_mesh.template | 7 ++++++ .../global_ocean/wave_mesh/scrip_file.py | 20 ++++++++++------- .../wave_mesh/scrip_file.template | 6 +++++ 6 files changed, 57 insertions(+), 27 deletions(-) create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.template create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.template create mode 100644 compass/ocean/tests/global_ocean/wave_mesh/scrip_file.template diff --git a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py index 1994b9054b..278e89a9fb 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py @@ -1,3 +1,6 @@ +from importlib.resources import read_text + +from jinja2 import Template from mpas_tools.logging import check_call from compass import Step @@ -26,16 +29,15 @@ def setup(self): super().setup() - f = open(f'{self.work_dir}/cull_waves_mesh.nml', 'w') - f.write('&inputs\n') - f.write(" waves_mesh_file = 'wave_base_mesh.nc'\n") - f.write(" ocean_mesh_file = 'ocean_culled_mesh.nc'\n") - f.write("/\n") - f.write("&output\n") - f.write(" waves_mesh_culled_vtk = 'wave_mesh_culled.vtk'\n") - f.write(" waves_mesh_culled_gmsh = 'wave_mesh_culled.msh'\n") - f.write("/\n") - f.close() + template = Template(read_text( + 'compass.ocean.tests.global_ocean.wave_mesh', + 'cull_mesh.template')) + + text = template.render() + text = f'{text}\n' + + with open(f'{self.work_dir}/cull_waves_mesh.nml', 'w') as file: + file.write(text) def run(self): diff --git a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.template b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.template new file mode 100644 index 0000000000..0c5dfe2b6b --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.template @@ -0,0 +1,8 @@ +&inputs + waves_mesh_file = 'wave_base_mesh.nc' + ocean_mesh_file = 'ocean_culled_mesh.nc' +/ +&output + waves_mesh_culled_vtk = 'wave_mesh_culled.vtk' + waves_mesh_culled_gmsh = 'wave_mesh_culled.msh' +/ diff --git a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py index e72d254227..483805377a 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.py @@ -1,3 +1,6 @@ +from importlib.resources import read_text + +from jinja2 import Template from mpas_tools.logging import check_call from compass import Step @@ -21,15 +24,15 @@ def setup(self): super().setup() - f = open(f'{self.work_dir}/rotate.nml', 'w') - f.write('&inputs\n') - f.write(" LON_POLE = -42.8906d0\n") - f.write(" LAT_POLE = 72.3200d0\n") - f.write(" wind_file = 'null'\n") - f.write(" mesh_file = 'wave_mesh_culled.msh'\n") - f.write(" mesh_file_out = 'wave_culled_mesh_RTD.msh'\n") - f.write("/\n") - f.close() + template = Template(read_text( + 'compass.ocean.tests.global_ocean.wave_mesh', + 'rotate_mesh.template')) + + text = template.render() + text = f'{text}\n' + + with open(f'{self.work_dir}/rotate.nml', 'w') as file: + file.write(text) def run(self): diff --git a/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.template b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.template new file mode 100644 index 0000000000..e7cd97d008 --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/rotate_mesh.template @@ -0,0 +1,7 @@ +&inputs + LON_POLE = -42.8906d0 + LAT_POLE = 72.3200d0 + wind_file = 'null' + mesh_file = 'wave_mesh_culled.msh' + mesh_file_out = 'wave_culled_mesh_RTD.msh' +/ diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index 0e60ccd98f..bb30f0f029 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -1,3 +1,6 @@ +from importlib.resources import read_text + +from jinja2 import Template from mpas_tools.logging import check_call from compass import Step @@ -21,14 +24,15 @@ def setup(self): super().setup() - f = open(f'{self.work_dir}/scrip.nml', 'w') - f.write('&inputs\n') - f.write(" waves_mesh_file = 'wave_mesh_culled.msh'\n") - f.write("/\n") - f.write('&outputs\n') - f.write(" waves_scrip_file = 'wave_mesh_scrip.nc'\n") - f.write("/\n") - f.close() + template = Template(read_text( + 'compass.ocean.tests.global_ocean.wave_mesh', + 'scrip_file.template')) + + text = template.render() + text = f'{text}\n' + + with open(f'{self.work_dir}/scrip.nml', 'w') as file: + file.write(text) def run(self): """ diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.template b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.template new file mode 100644 index 0000000000..c6d3582026 --- /dev/null +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.template @@ -0,0 +1,6 @@ +&inputs + waves_mesh_file = 'wave_mesh_culled.msh' +/ +&outputs + waves_scrip_file = 'wave_mesh_scrip.nc' +/ From c325b88d63a4b11021716ad61344a36cc6fde88a Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 14 Jun 2024 12:27:01 -0500 Subject: [PATCH 20/26] Create ocean scrip file and fix remapping step --- compass/ocean/tests/global_ocean/__init__.py | 2 +- .../tests/global_ocean/wave_mesh/__init__.py | 8 +-- .../global_ocean/wave_mesh/remap_files.py | 68 ++++++++++++------- .../global_ocean/wave_mesh/scrip_file.py | 11 ++- 4 files changed, 59 insertions(+), 30 deletions(-) diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index adbb797f41..9c7389c0f5 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -176,4 +176,4 @@ def _add_tests(self, mesh_names, high_res_topography=True, if include_wave_mesh: self.add_test_case( WaveMesh(test_group=self, ocean_mesh=mesh_test, - ocean_init=init_test, ocean_e3sm=e3sm_files_test)) + ocean_init=init_test)) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index d77b7a10e0..3a15bc8a2b 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -27,7 +27,7 @@ class WaveMesh(TestCase): Attributes ---------- """ - def __init__(self, test_group, ocean_mesh, ocean_init, ocean_e3sm): + def __init__(self, test_group, ocean_mesh, ocean_init): """ Create test case for creating a global MPAS-Ocean mesh @@ -54,12 +54,12 @@ def __init__(self, test_group, ocean_mesh, ocean_init, ocean_e3sm): self.add_step(cull_mesh_step) scrip_file_step = WavesScripFile(test_case=self, - wave_culled_mesh=cull_mesh_step) + wave_culled_mesh=cull_mesh_step, + ocean_mesh=ocean_init) self.add_step(scrip_file_step) remap_file_step = WavesRemapFiles(test_case=self, - wave_scrip=scrip_file_step, - ocean_e3sm=ocean_e3sm) + wave_scrip=scrip_file_step) self.add_step(remap_file_step) rotate_mesh_step = WavesRotateMesh(test_case=self, diff --git a/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py b/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py index 64afe1fbac..79f5e42a1d 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/remap_files.py @@ -1,4 +1,5 @@ from datetime import date +from distutils.spawn import find_executable from mpas_tools.logging import check_call @@ -9,7 +10,7 @@ class WavesRemapFiles(Step): """ A step for creating remapping files for wave mesh """ - def __init__(self, test_case, wave_scrip, ocean_e3sm, + def __init__(self, test_case, wave_scrip, name='remap_files', subdir=None): super().__init__(test_case=test_case, name=name, subdir=subdir) @@ -19,14 +20,15 @@ def __init__(self, test_case, wave_scrip, ocean_e3sm, filename='wave_scrip.nc', work_dir_target=f'{wave_scrip_file_path}/wave_mesh_scrip.nc') - ocean_scrip_file_path = ocean_e3sm.steps['scrip'].path + ocean_scrip_file_path = wave_scrip.path self.add_input_file( filename='ocean_scrip.nc', - work_dir_target=f'{ocean_scrip_file_path}/ocean_scrip.nc') + work_dir_target=f'{ocean_scrip_file_path}/ocean_mesh_scrip.nc') def make_remapping_files(self, grid1, grid2, grid1_shortname, grid2_shortname, datestamp, reg1, reg2, map_type, nprocs=1): + if map_type == 'conserve': map_abbrev = 'aave' elif map_type == 'bilinear': @@ -41,35 +43,53 @@ def make_remapping_files(self, grid1, grid2, print('map type not recognized') raise SystemExit(0) - flags = ' --ignore_unmapped --ignore_degenerate' - if reg1: - flags += ' --src_regional' - if reg2: - flags += ' --dst_regional' + exe = find_executable('ESMF_RegridWeightGen') + parallel_executable = self.config.get('parallel', + 'parallel_executable') + ntasks = self.ntasks + parallel_args = parallel_executable.split(' ') + if 'srun' in parallel_args: + parallel_args.extend(['-n', f'{ntasks}']) + else: # presume mpirun syntax + parallel_args.extend(['-np', f'{ntasks}']) + flags = ['--ignore_unmapped', '--ignore_degenerate'] + + # wave to ocean remap map_name = f'map_{grid1_shortname}_TO_{grid2_shortname}'\ f'_{map_abbrev}.{datestamp}.nc' - check_call(f'srun -n {nprocs}' - f' ESMF_RegridWeightGen --source {grid1}' - f' --destination {grid2}' - f' --method {map_type}' - f' --weight {map_name} {flags}', - logger=self.logger) - - flags = ' --ignore_unmapped --ignore_degenerate' + + args = [exe, + '--source', grid1, + '--destination', grid2, + '--method', map_type, + '--weight', map_name] + args.extend(flags) if reg1: - flags += ' --dst_regional' + args.append('--src_regional') if reg2: - flags += ' --src_regional' + args.append('--dst_regional') + cmd = parallel_args + args + check_call(cmd, logger=self.logger) + + # ocean to wave remap map_name = f'map_{grid2_shortname}_TO_{grid1_shortname}'\ f'_{map_abbrev}.{datestamp}.nc' - check_call(f'srun -n {nprocs}' - f' ESMF_RegridWeightGen --source {grid2}' - f' --destination {grid1}' - f' --method {map_type}' - f' --weight {map_name} {flags}', - logger=self.logger) + + args = [exe, + '--source', grid2, + '--destination', grid1, + '--method', map_type, + '--weight', map_name] + args.extend(flags) + if reg1: + args.append('--dst_regional') + if reg2: + args.append('--src_regional') + + cmd = parallel_args + args + check_call(cmd, logger=self.logger) def run(self): today = date.today() diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index bb30f0f029..83fd1a3438 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -2,6 +2,7 @@ from jinja2 import Template from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas from compass import Step @@ -10,11 +11,16 @@ class WavesScripFile(Step): """ A step for creating the scrip file for the wave mesh """ - def __init__(self, test_case, wave_culled_mesh, + def __init__(self, test_case, wave_culled_mesh, ocean_mesh, name='scrip_file', subdir=None): super().__init__(test_case=test_case, name=name, subdir=subdir) + ocean_mesh_path = ocean_mesh.steps['initial_state'].path + self.add_input_file( + filename='ocean_mesh.nc', + work_dir_target=f'{ocean_mesh_path}/initial_state.nc') + wave_culled_mesh_path = wave_culled_mesh.path self.add_input_file( filename='wave_mesh_culled.msh', @@ -38,4 +44,7 @@ def run(self): """ Create scrip files for wave mesh """ + local_filename = 'ocean_mesh_scrip.nc' + scrip_from_mpas('ocean_mesh.nc', local_filename) + check_call('ocean_scrip_wave_mesh', logger=self.logger) From 1fa223fe0f00594ee7f0bfd90f04a15415fda751 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Mon, 15 Jul 2024 11:49:24 -0500 Subject: [PATCH 21/26] Add UOST rotation --- .../tests/global_ocean/wave_mesh/__init__.py | 3 +- .../global_ocean/wave_mesh/uost_files.py | 206 +++++++++++++++++- .../global_ocean/wave_mesh/wave_mesh.cfg | 3 + 3 files changed, 204 insertions(+), 8 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 3a15bc8a2b..9fdc0f86c4 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/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): diff --git a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py b/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py index 76ddb5c357..7c1623650a 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/uost_files.py +++ b/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] diff --git a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg index 2ecf947586..12d8850d20 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg +++ b/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 From 849a837084fa5da639cdcf51b0d831d7e66a9883 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Tue, 23 Jul 2024 14:39:21 -0500 Subject: [PATCH 22/26] Add units to config file --- .../tests/global_ocean/wave_mesh/wave_mesh.cfg | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg index 12d8850d20..042a680489 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg +++ b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg @@ -1,13 +1,13 @@ [wave_mesh] -hfun_grid_spacing = 0.5 +hfun_grid_spacing = 0.5 # units of deg hfun_slope_lim = 0.15 -depth_threshold_refined = 1000.0 -distance_threshold_refined = 300.0 -depth_threshold_global = 1000.0 -distance_threshold_global = 300.0 -refined_res = 20000.0 -maxres = 225000.0 +depth_threshold_refined = 1000.0 # units of m +distance_threshold_refined = 300.0 # units of km +depth_threshold_global = 1000.0 # units of m +distance_threshold_global = 300.0 # units of km +refined_res = 20000.0 # units of m +maxres = 225000.0 # units of m ndir = 36 nfreq = 50 From 62a4aa5bd7159b345bcb9305215dedb1c9f715b8 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Thu, 5 Sep 2024 15:11:22 -0500 Subject: [PATCH 23/26] Change coastline data set from level 6 to level 5 - This means Antarctic coastline is between ocean and ice instead of ocean and grounding line --- compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index 726df345b0..6d96880e2e 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -195,7 +195,7 @@ def distance_to_shapefile_points(self, lon, lat, sphere_radius, reggrid=False): # Get coastline coordinates from shapefile - features = cfeature.GSHHSFeature(scale='l', levels=[1, 6]) + features = cfeature.GSHHSFeature(scale='l', levels=[1, 5]) pt_list = [] for feature in features.geometries(): From 1550d55ea67fd3885a9af4e54d186b08361fef06 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 6 Sep 2024 12:54:19 -0500 Subject: [PATCH 24/26] Incorporate ocean base mesh into wave base mesh creation - Also add the ability to (optionally) specify paths to existing base and culled ocean meshes in the config file --- .../tests/global_ocean/wave_mesh/__init__.py | 3 +- .../tests/global_ocean/wave_mesh/base_mesh.py | 105 +++++++++++++----- .../tests/global_ocean/wave_mesh/cull_mesh.py | 16 ++- .../global_ocean/wave_mesh/scrip_file.py | 16 ++- 4 files changed, 101 insertions(+), 39 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index 9fdc0f86c4..f75e67f7f6 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -45,7 +45,8 @@ def __init__(self, test_group, ocean_mesh, ocean_init): self.mesh_config_filename = 'wave_mesh.cfg' base_mesh_step = WavesBaseMesh(test_case=self, - ocean_mesh=ocean_init) + ocean_base_mesh=ocean_mesh, + ocean_culled_mesh=ocean_init) self.add_step(base_mesh_step) cull_mesh_step = WavesCullMesh(test_case=self, diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index 6d96880e2e..37c3a061b7 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt import numpy as np from mpas_tools.cime.constants import constants +from mpas_tools.mesh.cull import write_map_culled_to_base from netCDF4 import Dataset from scipy import interpolate, spatial @@ -16,21 +17,44 @@ class WavesBaseMesh(QuasiUniformSphericalMeshStep): """ A step for creating wave mesh based on an ocean mesh """ - def __init__(self, test_case, ocean_mesh, name='base_mesh', subdir=None): + def __init__(self, test_case, ocean_base_mesh, ocean_culled_mesh, + name='base_mesh', subdir=None): super().__init__(test_case=test_case, name=name, subdir=subdir, cell_width=None) - mesh_path = ocean_mesh.steps['initial_state'].path - self.add_input_file( - filename='ocean_mesh.nc', - work_dir_target=f'{mesh_path}/initial_state.nc') + self.ocean_base_mesh = ocean_base_mesh + self.ocean_culled_mesh = ocean_culled_mesh def setup(self): super().setup() self.opts.init_file = 'init.msh' + # Set links to base mesh + if self.config.has_option('wave_mesh', 'ocean_base_mesh'): + ocean_base_mesh_path = self.config.get('wave_mesh', + 'ocean_base_mesh') + else: + mesh_path = self.ocean_base_mesh.steps['base_mesh'].path + ocean_base_mesh_path = f'{mesh_path}/base_mesh.nc' + + self.add_input_file( + filename='ocean_base_mesh.nc', + work_dir_target=ocean_base_mesh_path) + + # Set links to culled mesh + if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): + ocean_culled_mesh_path = self.config.get('wave_mesh', + 'ocean_culled_mesh') + else: + mesh_path = self.ocean_culled_mesh.steps['initial_state'].path + ocean_culled_mesh_path = f'{mesh_path}/initial_state.nc' + + self.add_input_file( + filename='ocean_culled_mesh.nc', + work_dir_target=ocean_culled_mesh_path) + def build_cell_width_lat_lon(self): """ Create cell width array for this mesh on a regular latitude-longitude @@ -62,10 +86,13 @@ def build_cell_width_lat_lon(self): earth_radius = constants['SHR_CONST_REARTH'] / km cell_width = self.cell_widthVsLatLon(xlon, ylat, - earth_radius, 'ocean_mesh.nc') + earth_radius, + 'ocean_culled_mesh.nc') cell_width = cell_width / km - self.create_initial_points('ocean_mesh.nc', xlon, ylat, cell_width, + self.create_initial_points('ocean_base_mesh.nc', + 'ocean_culled_mesh.nc', + xlon, ylat, cell_width, earth_radius, self.opts.init_file) hfun_slope_lim = self.config.getfloat('wave_mesh', 'hfun_slope_lim') @@ -248,24 +275,29 @@ def distance_to_shapefile_points(self, lon, lat, return D - def create_initial_points(self, meshfile, lon, lat, hfunction, + def create_initial_points(self, base_mesh, culled_mesh, + lon, lat, hfunction, sphere_radius, outfile): - # Open MPAS mesh and get cell variables - nc_file = Dataset(meshfile, 'r') - lonCell = nc_file.variables['lonCell'][:] - latCell = nc_file.variables['latCell'][:] - nCells = lonCell.shape[0] + # Open MPAS culled mesh and get cell variables + nc_file_culled = Dataset(culled_mesh, 'r') + lonCell_culled = nc_file_culled.variables['lonCell'][:] + nCells_culled = lonCell_culled.shape[0] + + # Open MPAS base mesh and get cell variables + nc_file_base = Dataset(base_mesh, 'r') + lonCell_base = nc_file_base.variables['lonCell'][:] + latCell_base = nc_file_base.variables['latCell'][:] # Transform 0,360 range to -180,180 - idx, = np.where(lonCell > np.pi) - lonCell[idx] = lonCell[idx] - 2.0 * np.pi + idx, = np.where(lonCell_base > np.pi) + lonCell_base[idx] = lonCell_base[idx] - 2.0 * np.pi # Interpolate hfunction onto mesh cell centers hfun = interpolate.RegularGridInterpolator( (np.radians(lon), np.radians(lat)), hfunction.T) - mesh_pts = np.vstack((lonCell, latCell)).T + mesh_pts = np.vstack((lonCell_base, latCell_base)).T hfun_interp = hfun(mesh_pts) # Find cells in refined region of waves mesh @@ -277,34 +309,47 @@ def create_initial_points(self, meshfile, lon, lat, hfunction, # in the extra columns of cellsOnCell # in this case, these must be zeroed out # to correctly identify boundary cells - nEdgesOnCell = nc_file.variables['nEdgesOnCell'][:] - cellsOnCell = nc_file.variables['cellsOnCell'][:] - nz = np.zeros(cellsOnCell.shape, dtype=bool) - for i in range(nCells): - nz[i, 0:nEdgesOnCell[i]] = True - cellsOnCell[~nz] = 0 - nCellsOnCell = np.count_nonzero(cellsOnCell, axis=1) - is_boundary_cell = np.equal(nCellsOnCell, nEdgesOnCell) + nEdgesOnCell_culled = nc_file_culled.variables['nEdgesOnCell'][:] + cellsOnCell_culled = nc_file_culled.variables['cellsOnCell'][:] + nz = np.zeros(cellsOnCell_culled.shape, dtype=bool) + for i in range(nCells_culled): + nz[i, 0:nEdgesOnCell_culled[i]] = True + cellsOnCell_culled[~nz] = 0 + nCellsOnCell = np.count_nonzero(cellsOnCell_culled, axis=1) + is_boundary_cell = np.equal(nCellsOnCell, nEdgesOnCell_culled) idx_bnd, = np.where(is_boundary_cell == False) # noqa: E712 - # Force inclusion of all boundary cells - idx = np.union1d(idx, idx_bnd) + # Map from culled mesh cells to base mesh cells + write_map_culled_to_base(base_mesh, culled_mesh, + 'base_to_culled_map.nc') + map_file = Dataset('base_to_culled_map.nc', 'r') + culled_to_base = map_file.variables['mapCulledToBaseCell'][:] + base_idx_bnd = culled_to_base[idx_bnd] + + # Find cells in base mesh connected to culled mesh boundary cells + cellsOnCell_base = nc_file_base.variables['cellsOnCell'][:] + base_idx_xbnd = np.concatenate(cellsOnCell_base[base_idx_bnd]) + base_idx_xbnd = base_idx_xbnd - 1 + + # Force inclusion of all boundary cells, eliminating any duplicates + idx = np.union1d(idx, base_idx_bnd) + idx = np.union1d(idx, base_idx_xbnd) # Get coordinates of points - lon = lonCell[idx] - lat = latCell[idx] + lon = lonCell_base[idx] + lat = latCell_base[idx] + # Include north pole point lon = np.append(lon, 0.0) lat = np.append(lat, 0.5 * np.pi) - npt = lon.size - # Change to Cartesian coordinates x, y, z = self.lonlat2xyz(lon, lat, sphere_radius) # Get coordinates and ID into structured array # (for use with np.savetxt) pt_list = [] + npt = lon.size for i in range(npt): # ID of -1 specifies that node is fixed pt_list.append((x[i], y[i], z[i], -1)) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py index 278e89a9fb..f048a22c25 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py @@ -15,10 +15,7 @@ def __init__(self, test_case, ocean_mesh, wave_base_mesh, super().__init__(test_case=test_case, name=name, subdir=subdir) - culled_mesh_path = ocean_mesh.steps['initial_state'].path - self.add_input_file( - filename='ocean_culled_mesh.nc', - work_dir_target=f'{culled_mesh_path}/initial_state.nc') + self.ocean_mesh = ocean_mesh wave_base_mesh_path = wave_base_mesh.path self.add_input_file( @@ -29,6 +26,17 @@ def setup(self): super().setup() + if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): + culled_mesh_path = self.config.get('wave_mesh', + 'ocean_culled_mesh') + else: + mesh_path = self.ocean_mesh.steps['initial_state'].path + culled_mesh_path = f'{mesh_path}/initial_state.nc' + + self.add_input_file( + filename='ocean_culled_mesh.nc', + work_dir_target=culled_mesh_path) + template = Template(read_text( 'compass.ocean.tests.global_ocean.wave_mesh', 'cull_mesh.template')) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index 83fd1a3438..6c3a1e6509 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -16,10 +16,7 @@ def __init__(self, test_case, wave_culled_mesh, ocean_mesh, super().__init__(test_case=test_case, name=name, subdir=subdir) - ocean_mesh_path = ocean_mesh.steps['initial_state'].path - self.add_input_file( - filename='ocean_mesh.nc', - work_dir_target=f'{ocean_mesh_path}/initial_state.nc') + self.ocean_mesh = ocean_mesh wave_culled_mesh_path = wave_culled_mesh.path self.add_input_file( @@ -30,6 +27,17 @@ def setup(self): super().setup() + if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): + ocean_mesh_path = self.config.get('wave_mesh', + 'ocean_culled_mesh') + else: + mesh_path = self.ocean_mesh.steps['initial_state'].path + ocean_mesh_path = f'{mesh_path}/initial_state.nc' + + self.add_input_file( + filename='ocean_mesh.nc', + work_dir_target=ocean_mesh_path) + template = Template(read_text( 'compass.ocean.tests.global_ocean.wave_mesh', 'scrip_file.template')) From 66bb5c975ebbd370211e64fbe12b61c2f8955b07 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 6 Sep 2024 13:01:38 -0500 Subject: [PATCH 25/26] Clean up config file --- .../global_ocean/wave_mesh/wave_mesh.cfg | 45 +++++++++++++++---- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg index 042a680489..bec4b05128 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg +++ b/compass/ocean/tests/global_ocean/wave_mesh/wave_mesh.cfg @@ -1,13 +1,42 @@ [wave_mesh] -hfun_grid_spacing = 0.5 # units of deg +# Spacing of grid used to specify wave mesh resolution, units of deg +hfun_grid_spacing = 0.5 + +# Resolution gradient limit hfun_slope_lim = 0.15 -depth_threshold_refined = 1000.0 # units of m -distance_threshold_refined = 300.0 # units of km -depth_threshold_global = 1000.0 # units of m -distance_threshold_global = 300.0 # units of km -refined_res = 20000.0 # units of m -maxres = 225000.0 # units of m +# Ocean mesh resolution threshold for where depth_threshold_refined +# and distance_threshold_refined, criteria apply. This can be used +# to set more flexible threshold values for regionally refined ocean +# meshes, units of m. +refined_res = 20000.0 + +# Ocean resolution gets applied where depth is less than this threshold +# and where ocean resolution is less than refined_res, units of m +depth_threshold_refined = 1000.0 + +# Maximum distance from coast where ocean resolution gets applied +# where ocean resolutoin is less than refined_res, units of km +distance_threshold_refined = 300.0 + +# Ocean resolution gets applied where depth is less than this threshold, units of m +depth_threshold_global = 1000.0 + +# Maximum distance from ocast where ocean resolution gets applied, units of km +distance_threshold_global = 300.0 + +# Maximum resolution of waves mesh, units of m +maxres = 225000.0 + +# Number of spectral directions bins ndir = 36 -nfreq = 50 + +# Number of spectral frequency bins +nfreq = 36 + +# Path to existing base ocean mesh (optional) +# ocean_base_mesh = + +# Path to existing culled ocean mesh (optional) +# ocean_culled_mesh = From acc78b61a5db5a78c544cbc2f6e76accbee3328d Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 6 Sep 2024 16:15:03 -0500 Subject: [PATCH 26/26] Add standalone test for creating wave mesh from existing ocean mesh --- compass/ocean/tests/global_ocean/__init__.py | 3 ++- .../tests/global_ocean/wave_mesh/__init__.py | 8 +++++-- .../tests/global_ocean/wave_mesh/base_mesh.py | 24 ++++++++++++------- .../tests/global_ocean/wave_mesh/cull_mesh.py | 12 ++++++---- .../global_ocean/wave_mesh/scrip_file.py | 12 ++++++---- 5 files changed, 40 insertions(+), 19 deletions(-) diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index 9c7389c0f5..f07959f5fa 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -63,7 +63,8 @@ def __init__(self, mpas_core): # A test case for making E3SM support files from an existing mesh self.add_test_case(FilesForE3SM(test_group=self)) - # Create waves mesh + # Test cases for creating waves mesh + self.add_test_case(WaveMesh(test_group=self)) self._add_tests(mesh_names=['Icos'], include_wave_mesh=True) def _add_tests(self, mesh_names, high_res_topography=True, diff --git a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py index f75e67f7f6..037e08b99a 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/__init__.py @@ -27,7 +27,7 @@ class WaveMesh(TestCase): Attributes ---------- """ - def __init__(self, test_group, ocean_mesh, ocean_init): + def __init__(self, test_group, ocean_mesh=None, ocean_init=None): """ Create test case for creating a global MPAS-Ocean mesh @@ -38,7 +38,11 @@ def __init__(self, test_group, ocean_mesh, ocean_init): """ name = 'wave_mesh' - subdir = f'{ocean_mesh.mesh_name}/{name}' + if ocean_mesh is not None: + subdir = f'{ocean_mesh.mesh_name}/{name}' + else: + subdir = name + super().__init__(test_group=test_group, name=name, subdir=subdir) self.package = 'compass.ocean.tests.global_ocean.wave_mesh' diff --git a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py index 37c3a061b7..b6d9ec76e5 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/base_mesh.py @@ -32,24 +32,32 @@ def setup(self): self.opts.init_file = 'init.msh' # Set links to base mesh - if self.config.has_option('wave_mesh', 'ocean_base_mesh'): - ocean_base_mesh_path = self.config.get('wave_mesh', - 'ocean_base_mesh') - else: + if self.ocean_base_mesh is not None: mesh_path = self.ocean_base_mesh.steps['base_mesh'].path ocean_base_mesh_path = f'{mesh_path}/base_mesh.nc' + else: + if self.config.has_option('wave_mesh', 'ocean_base_mesh'): + ocean_base_mesh_path = self.config.get('wave_mesh', + 'ocean_base_mesh') + else: + raise ValueError('ocean_base_mesh option required ' + 'in wave_mesh section in cfg file') self.add_input_file( filename='ocean_base_mesh.nc', work_dir_target=ocean_base_mesh_path) # Set links to culled mesh - if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): - ocean_culled_mesh_path = self.config.get('wave_mesh', - 'ocean_culled_mesh') - else: + if self.ocean_culled_mesh is not None: mesh_path = self.ocean_culled_mesh.steps['initial_state'].path ocean_culled_mesh_path = f'{mesh_path}/initial_state.nc' + else: + if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): + ocean_culled_mesh_path = self.config.get('wave_mesh', + 'ocean_culled_mesh') + else: + raise ValueError('ocean_culled_mesh option required ' + 'in wave_mesh section in cfg file') self.add_input_file( filename='ocean_culled_mesh.nc', diff --git a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py index f048a22c25..4598be3681 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/cull_mesh.py @@ -26,12 +26,16 @@ def setup(self): super().setup() - if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): - culled_mesh_path = self.config.get('wave_mesh', - 'ocean_culled_mesh') - else: + if self.ocean_mesh is not None: mesh_path = self.ocean_mesh.steps['initial_state'].path culled_mesh_path = f'{mesh_path}/initial_state.nc' + else: + if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): + culled_mesh_path = self.config.get('wave_mesh', + 'ocean_culled_mesh') + else: + raise ValueError('ocean_culled_mesh option required ' + 'in wave_mesh section in cfg file') self.add_input_file( filename='ocean_culled_mesh.nc', diff --git a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py index 6c3a1e6509..73e23f88f7 100644 --- a/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py +++ b/compass/ocean/tests/global_ocean/wave_mesh/scrip_file.py @@ -27,12 +27,16 @@ def setup(self): super().setup() - if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): - ocean_mesh_path = self.config.get('wave_mesh', - 'ocean_culled_mesh') - else: + if self.ocean_mesh is not None: mesh_path = self.ocean_mesh.steps['initial_state'].path ocean_mesh_path = f'{mesh_path}/initial_state.nc' + else: + if self.config.has_option('wave_mesh', 'ocean_culled_mesh'): + ocean_mesh_path = self.config.get('wave_mesh', + 'ocean_culled_mesh') + else: + raise ValueError('ocean_culled_mesh option required ' + 'in wave_mesh section in cfg file') self.add_input_file( filename='ocean_mesh.nc',