Add subgrid cases to hurricane test group

compass/ocean/mesh/cull.py

@@ -14,7 +14,6 @@
 from mpas_tools.mesh.creation.sort_mesh import sort_mesh
 from mpas_tools.mesh.cull import cull_dataset, map_culled_to_base
 from mpas_tools.mesh.mask import compute_mpas_flood_fill_mask
-from mpas_tools.ocean import inject_bathymetry
 from mpas_tools.ocean.coastline_alteration import (
     add_critical_land_blockages,
     add_land_locked_cells_to_mask,
@@ -38,13 +37,9 @@ class CullMeshStep(Step):
     with_ice_shelf_cavities : bool
         Whether the mesh includes ice-shelf cavities
 
-    do_inject_bathymetry : bool
-        Whether to interpolate bathymetry from a data file so it
-        can be used as a culling criteria
-
     preserve_floodplain : bool
-        Whether to leave land cells in the mesh based on bathymetry
-        specified by do_inject_bathymetry
+        Whether to leave land cells in the mesh based on floodplain
+        specified by regionCellMasks
 
     unsmoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography
         A step for remapping topography. If provided, the remapped
@@ -57,7 +52,7 @@ class CullMeshStep(Step):
     """
 
     def __init__(self, test_case, base_mesh_step, with_ice_shelf_cavities,
-                 name='cull_mesh', subdir=None, do_inject_bathymetry=False,
+                 name='cull_mesh', subdir=None,
                  preserve_floodplain=False, unsmoothed_topo=None,
                  smoothed_topo=None):
         """
@@ -80,13 +75,10 @@ def __init__(self, test_case, base_mesh_step, with_ice_shelf_cavities,
         subdir : str, optional
             the subdirectory for the step.  The default is ``name``
 
-        do_inject_bathymetry : bool, optional
-            Whether to interpolate bathymetry from a data file so it
-            can be used as a culling criteria
 
         preserve_floodplain : bool, optional
-            Whether to leave land cells in the mesh based on bathymetry
-            specified by do_inject_bathymetry
+            Whether to leave land cells in the mesh based on floodplain
+            specified by regionCellMasks
 
         unsmoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography, optional
             A step for remapping topography. If provided, the remapped
@@ -110,7 +102,6 @@ def __init__(self, test_case, base_mesh_step, with_ice_shelf_cavities,
             self.add_output_file(filename='land_ice_mask.nc')
 
         self.with_ice_shelf_cavities = with_ice_shelf_cavities
-        self.do_inject_bathymetry = do_inject_bathymetry
         self.preserve_floodplain = preserve_floodplain
 
     def setup(self):
@@ -119,11 +110,6 @@ def setup(self):
         dependencies.
         """
         super().setup()
-        if self.do_inject_bathymetry:
-            self.add_input_file(
-                filename='earth_relief_15s.nc',
-                target='SRTM15_plus_earth_relief_15s.nc',
-                database='bathymetry_database')
 
         base_path = self.base_mesh_step.path
         base_filename = self.base_mesh_step.config.get(
@@ -180,7 +166,6 @@ def run(self):
         # only use progress bars if we're not writing to a log file
         use_progress_bar = self.log_filename is None
 
-        do_inject_bathymetry = self.do_inject_bathymetry
         preserve_floodplain = self.preserve_floodplain
 
         convert_to_cdf5 = config.getboolean('spherical_mesh',
@@ -200,9 +185,6 @@ def run(self):
                   latitude_threshold=latitude_threshold,
                   sweep_count=sweep_count)
 
-        if do_inject_bathymetry:
-            inject_bathymetry(mesh_file='culled_mesh.nc')
-
 
 def cull_mesh(with_cavities=False, with_critical_passages=False,
               custom_critical_passages=None, custom_land_blockages=None,
@@ -257,7 +239,7 @@ def cull_mesh(with_cavities=False, with_critical_passages=False,
         (``with_critical_passages=True``)
 
     preserve_floodplain : bool, optional
-        Whether to use the ``cellSeedMask`` field in the base mesh to preserve
+        Whether to use ``regionCellMasks`` field in the base mesh to preserve
         a floodplain at elevations above z=0
 
     logger : logging.Logger, optional
@@ -399,6 +381,7 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
         dsPreserve.append(dsCritPassMask)
 
     if preserve_floodplain:
+        # base mesh contains regionCellMasks field to indicate floodplain
         dsPreserve.append(dsBaseMesh)
 
     # fix land locked cells after adding critical land blockages, as these

compass/ocean/mesh/floodplain.py

@@ -1,5 +1,10 @@
-from mpas_tools.ocean import inject_bathymetry, inject_preserve_floodplain
+import netCDF4 as nc4
+import numpy as np
+from geometric_features import read_feature_collection
+from inpoly import inpoly2
+from mpas_tools.mesh.creation.signed_distance import _add_poly
 
+import compass.ocean.tests.tides.dem.dem_remap as dem_remap
 from compass.mesh.spherical import QuasiUniformSphericalMeshStep
 
 
@@ -44,6 +49,14 @@ def __init__(self, test_case, name='base_mesh', subdir=None,
                             target='SRTM15_plus_earth_relief_15s.nc',
                             database='bathymetry_database')
 
+        pixel_path = test_case.steps['pixel'].path
+        pixel_file = f'{pixel_path}/RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc'
+        self.add_input_file(
+            filename='bathy.nc',
+            work_dir_target=pixel_file)
+
+        self.bottomDepth_varname = 'bottomDepthObserved'
+
     def run(self):
         """
         Run this step of the test case
@@ -54,10 +67,119 @@ def run(self):
 
         mesh_filename = config.get('spherical_mesh', 'mpas_mesh_filename')
 
-        inject_bathymetry(mesh_file=mesh_filename)
+        dem_remap.dem_remap('bathy.nc', mesh_filename)
+
+        # Create new NetCDF variables in mesh file, if necessary
+        nc_mesh = nc4.Dataset(mesh_filename, 'r+')
+        nc_vars = nc_mesh.variables.keys()
+        if 'bottomDepthObserved' not in nc_vars:
+            nc_mesh.createVariable('bottomDepthObserved', 'f8', ('nCells'))
+
+        # Write to mesh file
+        nc_mesh.variables['bottomDepthObserved'][:] = \
+            nc_mesh.variables['bed_elevation'][:]
+        nc_mesh.close()
+
         if self.preserve_floodplain:
             floodplain_elevation = config.getfloat('spherical_mesh',
                                                    'floodplain_elevation')
-            inject_preserve_floodplain(
+            min_depth_outside_floodplain = config.getfloat(
+                'spherical_mesh',
+                'min_depth_outside_floodplain')
+            if config.has_option('spherical_mesh', 'floodplain_resolution'):
+                floodplain_resolution = config.getfloat(
+                    'spherical_mesh',
+                    'floodplain_resolution')
+            else:
+                floodplain_resolution = 1e10
+            if config.has_option('spherical_mesh', 'floodplain_geojson'):
+                floodplain_geojson = config.get(
+                    'spherical_mesh',
+                    'floodplain_geojson')
+                floodplain_region = read_feature_collection(floodplain_geojson)
+            else:
+                floodplain_region = None
+
+            self.inject_preserve_floodplain(
                 mesh_file=mesh_filename,
-                floodplain_elevation=floodplain_elevation)
+                floodplain_elevation=floodplain_elevation,
+                floodplain_resolution=floodplain_resolution,
+                floodplain_region=floodplain_region,
+                min_depth_outside_floodplain=min_depth_outside_floodplain)
+
+    def inject_preserve_floodplain(self, mesh_file, floodplain_elevation,
+                                   floodplain_resolution=1e10,
+                                   floodplain_region=None,
+                                   min_depth_outside_floodplain=5.0):
+
+        nc_mesh = nc4.Dataset(mesh_file, 'r+')
+        nc_vars = nc_mesh.variables.keys()
+
+        if 'regionCellMasks' not in nc_vars:
+            nc_mesh.createDimension('nRegions', 1)
+            nc_mesh.createVariable('regionCellMasks', 'i',
+                                   ('nCells', 'nRegions'))
+        nc_mesh.variables['regionCellMasks'][:] = 0.0
+
+        if 'transectCellMasks' not in nc_vars:
+            nc_mesh.createDimension('nTransects', 1)
+            nc_mesh.createVariable('transectCellMasks', 'i',
+                                   ('nCells', 'nTransects'))
+        nc_mesh.variables['transectCellMasks'][:] = 0.0
+
+        floodplain = self.find_floodplain(mesh_file, floodplain_elevation,
+                                          floodplain_resolution,
+                                          floodplain_region)
+
+        nc_mesh.variables['regionCellMasks'][:] = floodplain
+
+        nc_mesh.close()
+
+    def find_floodplain(self, mesh_file, floodplain_elevation,
+                        floodplain_resolution,
+                        floodplain_region):
+
+        nc_mesh = nc4.Dataset(mesh_file, 'r+')
+
+        bottomDepth = self.bottomDepth_varname
+        h = 2.0 * np.sqrt(nc_mesh.variables['areaCell'][:] / np.pi) / 1000.0
+        floodplain = np.logical_and(
+            h < floodplain_resolution,
+            nc_mesh.variables[bottomDepth][:] < floodplain_elevation)
+        if floodplain_region:
+            print("adding floodplain region")
+            nodes = list()
+            edges = list()
+            for feature in floodplain_region.features:
+                if feature['geometry']['type'] == 'Polygon':
+                    for poly in feature['geometry']['coordinates']:
+                        _add_poly(poly, edges, nodes)
+
+                elif feature['geometry']['type'] == 'MultiPolygon':
+                    for mpoly in feature['geometry']['coordinates']:
+                        for poly in mpoly:
+                            _add_poly(poly, edges, nodes)
+
+            nodes = np.array(nodes)
+            edges = np.array(edges)
+            print(nodes)
+            print(edges)
+
+            lon = np.degrees(nc_mesh.variables['lonCell'][:])
+            lon = np.mod(lon + 180.0, 360.0) - 180.0
+            lat = np.degrees(nc_mesh.variables['latCell'][:])
+
+            points = np.vstack([lon, lat]).T
+
+            floodplain_mask, _ = inpoly2(points, nodes, edges)
+            floodplain_mask = 1 * floodplain_mask
+            print(np.max(floodplain_mask))
+            print(np.min(floodplain_mask))
+
+        else:
+            floodplain_mask = 1
+        floodplain = floodplain_mask * floodplain
+
+        nc_mesh.close()
+
+        return floodplain

compass/ocean/tests/buttermilk_bay/buttermilk_bay.cfg

@@ -27,7 +27,8 @@ max_cells_per_core = 3000
 [buttermilk_bay_viz]
 
 # coordinates (in km) for timeseries plot
-points = [2.8, 0.53], [1.9, 1.66], [2.4, 3.029], [2.51, 3.027], [1.26, 1.56]
+#points = [2.8, 0.53], [1.9, 1.66], [2.4, 3.029], [2.51, 3.027], [1.26, 1.56]
+points = [2.8, 0.53], [1.9, 1.66], [3.0, 3.25], [2.6, 3.25], [1.26, 1.56]
 
 # generate contour plots at a specified interval between output timesnaps
 plot_interval = 1

compass/ocean/tests/buttermilk_bay/forward.py

@@ -1,5 +1,3 @@
-import time
-
 from compass.model import run_model
 from compass.step import Step
 
@@ -107,11 +105,27 @@ def get_dt(self):
         dt_per_m = config.getfloat('buttermilk_bay', 'dt_per_m')
 
         dt = dt_per_m * self.resolution
-        # https://stackoverflow.com/a/1384565/7728169
-        dt = time.strftime('%H:%M:%S', time.gmtime(dt))
+        dt = self.format_seconds(dt, 1)
 
         return dt
 
+    def format_seconds(self, seconds, decimal_places=2):
+        # Handle negative values
+        is_negative = seconds < 0
+        seconds = abs(seconds)
+
+        # Extract hours, minutes, and seconds
+        hours = int(seconds // 3600)
+        remaining = seconds % 3600
+        minutes = int(remaining // 60)
+        secs = remaining % 60
+
+        # Format the string with specified decimal places
+        sign = '-' if is_negative else ''
+        time = (f"{sign}{hours:02d}:{minutes:02d}:"
+                f"{secs:0{3 + decimal_places}.{decimal_places}f}")
+        return time
+
     def _get_resources(self):
         """ get the these properties from the config options """
         config = self.config

compass/ocean/tests/buttermilk_bay/initial_state.py

@@ -1,3 +1,4 @@
+import numpy as np
 from mpas_tools.io import write_netcdf
 from mpas_tools.mesh.conversion import convert, cull
 from mpas_tools.planar_hex import make_planar_hex_mesh
@@ -23,6 +24,7 @@ def __init__(self, test_case, name, resolution,
         """
         self.coord_type = coord_type
         self.resolution = resolution
+        self.wetdry = wetdry
 
         super().__init__(test_case=test_case, name=name, ntasks=1,
                          min_tasks=1, openmp_threads=1)
@@ -40,14 +42,19 @@ def __init__(self, test_case, name, resolution,
         self.add_input_file(
             filename='buttermilk_bathy.nc',
             target='buttermilk_bathy.nc',
-            database='bathymetry_database')
+            database='buttermilk_bay')
 
         for file in ['base_mesh.nc', 'culled_mesh.nc', 'culled_graph.info',
                      'ocean.nc', 'init_mode_forcing_data.nc']:
             self.add_output_file(file)
 
         self.add_model_as_input()
 
+    def setup(self):
+        level = int(self.resolution * np.tan(np.pi / 6.0)) + 1
+        option = {'config_subgrid_refinement_level': f'{level}'}
+        self.add_namelist_options(option, mode='init')
+
     def run(self):
         """
         Run this step of the test case
@@ -62,6 +69,7 @@ def run(self):
         else:
             vert_levels = config.getint('vertical_grid', 'vert_levels')
             options = {'config_buttermilk_bay_vert_levels': f'{vert_levels}'}
+
         self.update_namelist_at_runtime(options)
 
         # Determine mesh parameters

compass/ocean/tests/buttermilk_bay/namelist.forward

@@ -17,6 +17,8 @@ config_mom_del2 = 2.0e1
 config_use_tidal_forcing = .true.
 config_tidal_forcing_type = 'direct'
 config_tidal_forcing_model = 'monochromatic'
+config_use_variable_drag = .true.
+config_implicit_bottom_drag_type = 'mannings'
 config_zero_drying_velocity_ramp = .true.
 config_zero_drying_velocity_ramp_hmin = 0.05
 config_zero_drying_velocity_ramp_hmax = 0.10

compass/ocean/tests/buttermilk_bay/namelist.init

@@ -4,6 +4,15 @@ config_use_wetting_drying = .true.
 config_drying_min_cell_height = 0.05
 config_write_cull_cell_mask = .false.
 config_Buttermilk_bay_vert_levels = 1
-config_Buttermilk_bay_topography_file = 'buttermilk_bathy.nc'
 config_write_cull_cell_mask = .true.
 config_use_tidal_forcing = .true.
+config_use_variable_drag = .true.
+config_subgrid_topography_file = 'buttermilk_bathy.nc'
+config_use_subgrid_wetting_drying = .true.
+config_subgrid_table_levels = -1
+config_Buttermilk_bay_subgrid_table_levels = 100
+config_subgrid_ssh_table_min = -3.0
+config_subgrid_ssh_table_max = 3.0
+config_subgrid_refinement_level = 10
+config_subgrid_use_thin_layer = .false.
+config_subgrid_edge_bathymetry_max_pixel = .true.

compass/ocean/tests/buttermilk_bay/namelist.subgrid.init

@@ -1,6 +1,6 @@
 config_use_subgrid_wetting_drying = .true.
 config_subgrid_table_levels = -1
 config_Buttermilk_bay_subgrid_table_levels = 100
-config_Buttermilk_bay_subgrid_refinement_level = 10
-config_Buttermilk_bay_subgrid_use_thin_layer = .false.
-config_Buttermilk_bay_subgrid_edge_bathymetry_max_pixel = .true.
+config_subgrid_refinement_level = 10
+config_subgrid_use_thin_layer = .false.
+config_subgrid_edge_bathymetry_max_pixel = .true.

compass/ocean/tests/buttermilk_bay/streams.forward

@@ -48,4 +48,12 @@
     <var name="tidalInputMask"/>
 </stream>
 
+<stream name="bottomDrag"
+        type="input"
+        filename_template="init.nc"
+        input_interval="initial_only">
+
+        <var name="bottomDrag"/>
+</stream>
+
 </streams>