Merge pull request #716 from cbegeman/ocn-enhance-drying-slope-init

Replace init mode in drying slope test cases

compass/ocean/tests/drying_slope/analysis.py

@@ -91,7 +91,10 @@ def _compute_rmse(self, ds, t):
         ds = ds.drop_vars(np.setdiff1d([j for j in ds.variables],
                                        ['yCell', 'ssh']))
         ds = ds.isel(Time=tidx)
-        x_mpas = ds.yCell.values / 1000.0
+        drying_length = self.config.getfloat('drying_slope', 'ly_analysis')
+        drying_length = drying_length * 1e3
+        x_offset = np.max(ds.yCell.values) - drying_length
+        x_mpas = (ds.yCell.values - x_offset) / 1000.0
         ssh_mpas = ds.ssh.values
         # Interpolate mpas solution to the points at which we have an exact
         # solution

compass/ocean/tests/drying_slope/convergence/__init__.py

@@ -92,8 +92,7 @@ def _setup_steps(self, config, subdir, method):
             init_name = f'initial_state_{res_name}'
             self.add_step(InitialState(test_case=self,
                                        name=init_name,
-                                       resolution=resolution,
-                                       coord_type=self.coord_type))
+                                       resolution=resolution))
             ntasks = max(min_tasks,
                          int(ceil(ntasks_baseline / resolution**2.)))
             forward_step = Forward(test_case=self, resolution=resolution,
@@ -124,3 +123,9 @@ def validate(self):
                 res_name = f'{int(resolution)}km'
             compare_variables(test_case=self, variables=variables,
                               filename1=f'forward_{res_name}/output.nc')
+
+    def configure(self):
+        """
+        Change config options as needed
+        """
+        self.config.set('vertical_grid', 'coord_type', self.coord_type)

compass/ocean/tests/drying_slope/decomp/__init__.py

@@ -46,8 +46,7 @@ def __init__(self, test_group, resolution, coord_type, method):
         subdir = f'{coord_type}/{method}/{res_name}/{name}'
         super().__init__(test_group=test_group, name=name,
                          subdir=subdir)
-        self.add_step(InitialState(test_case=self, resolution=resolution,
-                                   coord_type=coord_type))
+        self.add_step(InitialState(test_case=self, resolution=resolution))
 
         if coord_type == 'single_layer':
             damping_coeff = None

compass/ocean/tests/drying_slope/default/__init__.py

@@ -50,8 +50,7 @@ def __init__(self, test_group, resolution, coord_type, method):
         subdir = f'{coord_type}/{method}/{res_name}/{name}'
         super().__init__(test_group=test_group, name=name,
                          subdir=subdir)
-        self.add_step(InitialState(test_case=self, resolution=resolution,
-                                   coord_type=coord_type))
+        self.add_step(InitialState(test_case=self, resolution=resolution))
         damping_coeffs = None
         config = CompassConfigParser()
         config.add_from_package('compass.ocean.tests.drying_slope',
@@ -99,3 +98,9 @@ def validate(self):
         else:
             compare_variables(test_case=self, variables=variables,
                               filename1='forward/output.nc')
+
+    def configure(self):
+        """
+        Change config options as needed
+        """
+        self.config.set('vertical_grid', 'coord_type', self.coord_type)

compass/ocean/tests/drying_slope/drying_slope.cfg

@@ -1,18 +1,57 @@
 # Options related to the vertical grid
 [vertical_grid]
 
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+bottom_depth = 10.
+
 # Number of vertical levels
 vert_levels = 10
 
-# Thickness of each layer in the thin film region
+# Thickness of each layer in the thin film region in m
 thin_film_thickness = 1.0e-3
 
+# Whether to use "partial" or "full", or "None" to not alter the topography
+partial_cell_type = None
+
+# The minimum fraction of a layer for partial cells
+min_pc_fraction = 0.1
+
 # config options for drying slope test cases
 [drying_slope]
 
 # the number of grid cells in x
 nx = 6
 
+# Length over which wetting and drying actually occur in km
+ly_analysis = 25.
+
+# Domain length in km
+ly = 30.
+
+# Bottom depth at the right side of the domain in m
+right_bottom_depth = 10.
+
+# Bottom depth at the left side of the domain in m
+left_bottom_depth = 0.
+
+# Plug width as a fraction of the domain
+plug_width_frac = 0.0
+
+# Plug temperature
+plug_temperature = 20.0
+
+# Background temperature
+background_temperature = 20.0
+
+# Background salinity
+background_salinity = 35.0
+
+# Coriolis parameter
+coriolis_parameter = 0.0
+
 # time step in s per km of horizontal resolution
 dt_per_km = 30
 
@@ -25,6 +64,7 @@ min_tasks = 1
 # config options for visualizing drying slope ouptut
 [drying_slope_convergence]
 
+# horizontal resolutions in km
 resolutions = 0.25, 0.5, 1, 2
 
 # config options for visualizing drying slope ouptut

compass/ocean/tests/drying_slope/forward.py

@@ -74,10 +74,10 @@ def __init__(self, test_case, resolution, name='forward', subdir=None,
                             target=f'{input_path}/culled_mesh.nc')
 
         self.add_input_file(filename='init.nc',
-                            target=f'{input_path}/ocean.nc')
+                            target=f'{input_path}/initial_state.nc')
 
         self.add_input_file(filename='forcing.nc',
-                            target=f'{input_path}/init_mode_forcing_data.nc')
+                            target=f'{input_path}/forcing.nc')
 
         self.add_input_file(filename='graph.info',
                             target=f'{input_path}/culled_graph.info')

compass/ocean/tests/drying_slope/initial_state.py

@@ -1,8 +1,9 @@
+import xarray as xr
 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
 
-from compass.model import run_model
+from compass.ocean.vertical import init_vertical_coord
 from compass.step import Step
 
 
@@ -12,7 +13,7 @@ class InitialState(Step):
     cases
     """
     def __init__(self, test_case, resolution, name='initial_state',
-                 coord_type='sigma'):
+                 baroclinic=False):
         """
         Create the step
 
@@ -24,18 +25,10 @@ def __init__(self, test_case, resolution, name='initial_state',
         super().__init__(test_case=test_case, name=name, ntasks=1,
                          min_tasks=1, openmp_threads=1)
 
-        self.coord_type = coord_type
-
         self.resolution = resolution
 
-        self.add_namelist_file('compass.ocean.tests.drying_slope',
-                               'namelist.init', mode='init')
-
-        self.add_streams_file('compass.ocean.tests.drying_slope',
-                              'streams.init', mode='init')
-
         for file in ['base_mesh.nc', 'culled_mesh.nc', 'culled_graph.info',
-                     'ocean.nc']:
+                     'initial_state.nc', 'forcing.nc']:
             self.add_output_file(file)
 
         self.add_model_as_input()
@@ -46,42 +39,96 @@ def run(self):
         """
         config = self.config
         logger = self.logger
+        resolution = self.resolution
 
+        # Fetch config options
         section = config['vertical_grid']
-        coord_type = self.coord_type
-        thin_film_thickness = section.getfloat('thin_film_thickness') + 1.0e-9
-        if coord_type == 'single_layer':
-            options = {'config_tidal_boundary_vert_levels': '1',
-                       'config_drying_min_cell_height':
-                       f'{thin_film_thickness}'}
-            self.update_namelist_at_runtime(options)
-        else:
-            vert_levels = section.getint('vert_levels')
-            options = {'config_tidal_boundary_vert_levels': f'{vert_levels}',
-                       'config_tidal_boundary_layer_type': f"'{coord_type}'",
-                       'config_drying_min_cell_height':
-                       f'{thin_film_thickness}'}
-            self.update_namelist_at_runtime(options)
+        thin_film_thickness = section.getfloat('thin_film_thickness')
+        vert_levels = section.getint('vert_levels')
+
+        section = config['drying_slope']
+        nx = section.getint('nx')
+        domain_length = section.getfloat('ly') * 1e3
+        drying_length = section.getfloat('ly_analysis') * 1e3
+        plug_width_frac = section.getfloat('plug_width_frac')
+        right_bottom_depth = section.getfloat('right_bottom_depth')
+        left_bottom_depth = section.getfloat('left_bottom_depth')
+        plug_temperature = section.getfloat('plug_temperature')
+        background_temperature = section.getfloat('background_temperature')
+        background_salinity = section.getfloat('background_salinity')
+        coriolis_parameter = section.getfloat('coriolis_parameter')
+
+        # Check config options
+        if domain_length < drying_length:
+            raise ValueError('Domain is not long enough to capture wetting '
+                             'front')
+        if right_bottom_depth < left_bottom_depth:
+            raise ValueError('Right boundary must be deeper than left '
+                             'boundary')
 
         # Determine mesh parameters
-        nx = config.getint('drying_slope', 'nx')
-        ny = round(28 / self.resolution)
-        if self.resolution < 1.:
+        dc = 1e3 * resolution
+        ny = round(domain_length / dc)
+        # This is just for consistency with previous implementations and could
+        # be removed
+        if resolution < 1.:
             ny += 2
-        dc = 1e3 * self.resolution
+        ny = 2 * round(ny / 2)
 
         logger.info(' * Make planar hex mesh')
-        dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=False,
-                                      nonperiodic_y=True)
+        ds_mesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc,
+                                       nonperiodic_x=False,
+                                       nonperiodic_y=True)
         logger.info(' * Completed Make planar hex mesh')
-        write_netcdf(dsMesh, 'base_mesh.nc')
+        write_netcdf(ds_mesh, 'base_mesh.nc')
 
         logger.info(' * Cull mesh')
-        dsMesh = cull(dsMesh, logger=logger)
+        ds_mesh = cull(ds_mesh, logger=logger)
         logger.info(' * Convert mesh')
-        dsMesh = convert(dsMesh, graphInfoFileName='culled_graph.info',
-                         logger=logger)
+        ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info',
+                          logger=logger)
         logger.info(' * Completed Convert mesh')
-        write_netcdf(dsMesh, 'culled_mesh.nc')
-        run_model(self, namelist='namelist.ocean',
-                  streams='streams.ocean')
+        write_netcdf(ds_mesh, 'culled_mesh.nc')
+
+        ds = ds_mesh.copy()
+        ds_forcing = ds_mesh.copy()
+
+        y_min = ds_mesh.yCell.min()
+        y_max = ds_mesh.yCell.max()
+        dc_edge_min = ds_mesh.dcEdge.min()
+
+        y_cell = ds.yCell
+        max_level_cell = vert_levels
+        bottom_depth = (right_bottom_depth - (y_max - y_cell) / drying_length *
+                        (right_bottom_depth - left_bottom_depth))
+        ds['bottomDepth'] = bottom_depth
+        # Set the water column to dry everywhere
+        ds['ssh'] = -bottom_depth + thin_film_thickness * max_level_cell
+        # We don't use config_tidal_forcing_monochromatic_baseline because the
+        # default value doesn't alter the initial state
+        init_vertical_coord(config, ds)
+
+        plug_width = domain_length * plug_width_frac
+        y_plug_boundary = y_min + plug_width
+        ds['temperature'] = xr.where(y_cell < y_plug_boundary,
+                                     plug_temperature, background_temperature)
+        ds['tracer1'] = xr.where(y_cell < y_plug_boundary, 1.0, 0.0)
+        ds['salinity'] = background_salinity * xr.ones_like(y_cell)
+        normalVelocity = xr.zeros_like(ds_mesh.xEdge)
+        normalVelocity, _ = xr.broadcast(normalVelocity, ds.refBottomDepth)
+        normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels')
+        normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0)
+        ds['normalVelocity'] = normalVelocity
+        ds['fCell'] = coriolis_parameter * xr.ones_like(ds.xCell)
+        ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge)
+        ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex)
+
+        write_netcdf(ds, 'initial_state.nc')
+
+        # Define the tidal boundary condition over 1-cell width
+        y_tidal_boundary = y_max - dc_edge_min / 2.
+        tidal_forcing_mask = xr.where(y_cell > y_tidal_boundary, 1.0, 0.0)
+        if tidal_forcing_mask.sum() <= 0:
+            raise ValueError('Input mask for tidal case is not set!')
+        ds_forcing['tidalInputMask'] = tidal_forcing_mask
+        write_netcdf(ds_forcing, 'forcing.nc')

compass/ocean/tests/drying_slope/loglaw/__init__.py

@@ -47,8 +47,7 @@ def __init__(self, test_group, resolution, coord_type, method):
         subdir = f'{coord_type}/{method}/{res_name}/{name}'
         super().__init__(test_group=test_group, name=name,
                          subdir=subdir)
-        self.add_step(InitialState(test_case=self, coord_type=coord_type,
-                                   resolution=resolution))
+        self.add_step(InitialState(test_case=self, resolution=resolution))
         config = CompassConfigParser()
         config.add_from_package('compass.ocean.tests.drying_slope',
                                 'drying_slope.cfg')
@@ -72,3 +71,9 @@ def validate(self):
         variables = ['layerThickness', 'normalVelocity']
         compare_variables(test_case=self, variables=variables,
                           filename1='forward/output.nc')
+
+    def configure(self):
+        """
+        Change config options as needed
+        """
+        self.config.set('vertical_grid', 'coord_type', self.coord_type)

compass/ocean/tests/drying_slope/namelist.forward

@@ -20,3 +20,4 @@ config_thickness_flux_type='upwind'
 config_use_cvmix=.false.
 config_use_debugTracers=.false.
 config_check_ssh_consistency=.true.
+config_disable_tr_all_tend=.true.