Add a step for smoothing topograpy

compass/ocean/mesh/cull.py

@@ -46,15 +46,20 @@ class CullMeshStep(Step):
         Whether to leave land cells in the mesh based on bathymetry
         specified by do_inject_bathymetry
 
-    remap_topography : compass.ocean.mesh.remap_topography.RemapTopography
+    unsmoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography
         A step for remapping topography. If provided, the remapped
         topography is used to determine the land mask
 
+    smoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography
+        A step for remapping topography. If provided, the remapped
+        topography is culled for subsequent use in ocean initial conditions
+
     """
 
     def __init__(self, test_case, base_mesh_step, with_ice_shelf_cavities,
                  name='cull_mesh', subdir=None, do_inject_bathymetry=False,
-                 preserve_floodplain=False, remap_topography=None):
+                 preserve_floodplain=False, unsmoothed_topo=None,
+                 smoothed_topo=None):
         """
         Create a new step
 
@@ -83,14 +88,19 @@ def __init__(self, test_case, base_mesh_step, with_ice_shelf_cavities,
             Whether to leave land cells in the mesh based on bathymetry
             specified by do_inject_bathymetry
 
-        remap_topography : compass.ocean.mesh.remap_topography.RemapTopography, optional
+        unsmoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography, optional
             A step for remapping topography. If provided, the remapped
             topography is used to determine the land mask
+
+        smoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography, optional
+            A step for remapping topography. If provided, the remapped
+            topography is culled for subsequent use in ocean initial conditions
         """  # noqa: E501
         super().__init__(test_case, name=name, subdir=subdir,
                          cpus_per_task=None, min_cpus_per_task=None)
         self.base_mesh_step = base_mesh_step
-        self.remap_topography = remap_topography
+        self.unsmoothed_topo = unsmoothed_topo
+        self.smoothed_topo = smoothed_topo
 
         for file in ['culled_mesh.nc', 'culled_graph.info',
                      'critical_passages_mask_final.nc']:
@@ -121,10 +131,16 @@ def setup(self):
         target = os.path.join(base_path, base_filename)
         self.add_input_file(filename='base_mesh.nc', work_dir_target=target)
 
-        if self.remap_topography is not None:
-            topo_path = self.remap_topography.path
+        if self.unsmoothed_topo is not None:
+            if self.smoothed_topo is None:
+                self.smoothed_topo = self.unsmoothed_topo
+            topo_path = self.unsmoothed_topo.path
             target = os.path.join(topo_path, 'topography_remapped.nc')
-            self.add_input_file(filename='topography.nc',
+            self.add_input_file(filename='unsmoothed_topography.nc',
+                                work_dir_target=target)
+            topo_path = self.smoothed_topo.path
+            target = os.path.join(topo_path, 'topography_remapped.nc')
+            self.add_input_file(filename='smoothed_topography.nc',
                                 work_dir_target=target)
             self.add_output_file('topography_culled.nc')
             self.add_output_file('map_culled_to_base.nc')
@@ -290,14 +306,14 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
     netcdf_format = mpas_tools.io.default_format
     netcdf_engine = mpas_tools.io.default_engine
 
-    has_remapped_topo = os.path.exists('topography.nc')
+    has_remapped_topo = os.path.exists('unsmoothed_topography.nc')
     if with_cavities and not has_remapped_topo:
         raise ValueError('Mesh culling with caviites must be from '
                          'remapped topography.')
 
     if has_remapped_topo:
         _land_mask_from_topo(with_cavities,
-                             topo_filename='topography.nc',
+                             topo_filename='unsmoothed_topography.nc',
                              mask_filename='land_mask.nc')
     else:
         _land_mask_from_geojson(with_cavities=with_cavities,
@@ -478,7 +494,7 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
 def _cull_topo(with_cavities, process_count, logger, latitude_threshold,
                sweep_count, ds_preserve):
 
-    ds_topo = xr.open_dataset('topography.nc')
+    ds_topo = xr.open_dataset('smoothed_topography.nc')
     ds_base = xr.open_dataset('base_mesh.nc')
 
     ds_culled = xr.open_dataset('culled_mesh.nc')
@@ -488,6 +504,15 @@ def _cull_topo(with_cavities, process_count, logger, latitude_threshold,
     write_netcdf(ds_map_culled_to_base, 'map_culled_to_base.nc')
 
     if with_cavities:
+        ds_unsmoothed_topo = xr.open_dataset('unsmoothed_topography.nc')
+
+        # use fractions from unsmoothed topography so we don't have
+        # different masks for different amounts of smoothing.  This
+        # allows us to use the same E3SM mapping files for different smoothing
+        for var in ['landIceFracObserved', 'landIceGroundedFracObserved',
+                    'landIceFloatingFracObserved']:
+            ds_topo[var] = ds_unsmoothed_topo[var]
+
         _flood_fill_and_add_land_ice_mask(ds_topo, ds_base,
                                           ds_map_culled_to_base, ds_preserve,
                                           logger, latitude_threshold,

compass/ocean/mesh/remap_topography.cfg

@@ -14,7 +14,7 @@ description = Bathymetry is from GEBCO 2023, combined with BedMachine
               Antarctica v3 around Antarctica.
 
 # the target and minimum number of MPI tasks to use in remapping
-ntasks = 1280
+ntasks = 640
 min_tasks = 256
 
 # remapping method {'bilinear', 'neareststod', 'conserve'}
@@ -29,5 +29,5 @@ renorm_threshold = 0.01
 ice_density = 910.0
 
 # smoothing parameters (no smoothing by default)
-expand_dist = 0
+expand_distance = 0
 expand_factor = 1

compass/ocean/mesh/remap_topography.py

@@ -8,6 +8,7 @@
 from mpas_tools.logging import check_call
 from pyremap import MpasCellMeshDescriptor
 
+from compass.io import symlink
 from compass.parallel import run_command
 from compass.step import Step
 
@@ -24,11 +25,17 @@ class RemapTopography(Step):
 
     mesh_name : str
         The name of the MPAS mesh to include in the mapping file
+
+    smoothing : bool, optional
+        Whether smoothing will be applied as part of the remapping
+
+    unsmoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography
+        A step with unsmoothed topography
     """
 
     def __init__(
         self, test_case, base_mesh_step, name='remap_topography', subdir=None,
-        mesh_name='MPAS_mesh',
+        mesh_name='MPAS_mesh', smoothing=False, unsmoothed_topo=None
     ):
         """
         Create a new step
@@ -49,11 +56,19 @@ def __init__(
 
         mesh_name : str, optional
             The name of the MPAS mesh to include in the mapping file
-        """
+
+        smoothing : bool, optional
+            Whether smoothing will be applied as part of the remapping
+
+        unsmoothed_topo : compass.ocean.mesh.remap_topography.RemapTopography, optional
+            A step with unsmoothed topography
+        """  # noqa: E501
         super().__init__(test_case, name=name, subdir=subdir,
                          ntasks=None, min_tasks=None)
         self.base_mesh_step = base_mesh_step
         self.mesh_name = mesh_name
+        self.smoothing = smoothing
+        self.unsmoothed_topo = unsmoothed_topo
 
         self.add_output_file(filename='topography_remapped.nc')
 
@@ -108,6 +123,11 @@ def run(self):
         """
         Run this step of the test case
         """
+        super().run()
+        if self._symlink_unsmoothed():
+            # we symlinked to the unsmoothed topography and we're done!
+            return
+
         config = self.config
         weight_generator = config.get('remap_topography', 'weight_generator')
 
@@ -124,25 +144,56 @@ def run(self):
         self._remap_to_target()
         self._modify_remapped_bathymetry()
 
+    def _symlink_unsmoothed(self):
+        """
+        If we are smoothing but no smoothing was actually requested, symlink
+        to the unsmoothed topography
+        """
+        if not self.smoothing or self.unsmoothed_topo is None:
+            # there's no unsmoothed topogrpahy yet
+            return False
+
+        config = self.config
+        section = config['remap_topography']
+        expand_distance = section.getfloat('expand_distance')
+        expand_factor = section.getfloat('expand_factor')
+
+        if expand_distance != 0. or expand_factor != 1.:
+            # we're doing some smoothing!
+            return False
+
+        # we already have unsmoothed topography and we're not doing
+        # smoothing so we can just symlink the unsmoothed results
+        out_filename = 'topography_remapped.nc'
+        unsmoothed_path = self.unsmoothed_topo.work_dir
+        target = os.path.join(unsmoothed_path, out_filename)
+        symlink(target, out_filename)
+
+        return True
+
     def _create_target_scrip_file(self):
         """
         Create target SCRIP file from MPAS mesh file.
         """
         logger = self.logger
         logger.info('Create source SCRIP file')
 
-        config = self.config
-        section = config['remap_topography']
-        expand_dist = section.getfloat('expand_dist')
-        expand_factor = section.getfloat('expand_factor')
+        if self.smoothing:
+            config = self.config
+            section = config['remap_topography']
+            expand_distance = section.getfloat('expand_distance')
+            expand_factor = section.getfloat('expand_factor')
+        else:
+            expand_distance = 0.
+            expand_factor = 1.
 
         descriptor = MpasCellMeshDescriptor(
             filename='base_mesh.nc',
             mesh_name=self.mesh_name,
         )
         descriptor.to_scrip(
             'target.scrip.nc',
-            expand_dist=expand_dist,
+            expand_dist=expand_distance,
             expand_factor=expand_factor,
         )
 

compass/ocean/suites/pr_to_cache.txt

@@ -0,0 +1,13 @@
+ocean/global_convergence/qu/cosine_bell
+
+ocean/global_ocean/IcoswISC240/mesh
+ocean/global_ocean/IcoswISC240/WOA23/init
+ocean/global_ocean/IcoswISC240/WOA23/performance_test
+
+ocean/global_ocean/Icos/mesh
+ocean/global_ocean/Icos/WOA23/init
+ocean/global_ocean/Icos/WOA23/performance_test
+
+ocean/global_ocean/IcoswISC/mesh
+ocean/global_ocean/IcoswISC/WOA23/init
+ocean/global_ocean/IcoswISC/WOA23/performance_test

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

@@ -148,21 +148,53 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
         self.add_step(base_mesh_step)
 
         if mali_ais_topo is None:
-            remap_step = RemapTopography(test_case=self,
-                                         base_mesh_step=base_mesh_step,
-                                         mesh_name=mesh_name)
+            unsmoothed_topo = RemapTopography(
+                test_case=self,
+                base_mesh_step=base_mesh_step,
+                name='remap_topo_unsmoothed',
+                mesh_name=mesh_name,
+                smoothing=False)
+
+            self.add_step(unsmoothed_topo)
+
+            smoothed_topo = RemapTopography(
+                test_case=self,
+                base_mesh_step=base_mesh_step,
+                name='remap_topo_smoothed',
+                mesh_name=mesh_name,
+                smoothing=True,
+                unsmoothed_topo=unsmoothed_topo)
+
+            self.add_step(smoothed_topo)
+
         else:
-            remap_step = RemapMaliTopography(
-                test_case=self, base_mesh_step=base_mesh_step,
-                mesh_name=mesh_name, mali_ais_topo=mali_ais_topo,
+            unsmoothed_topo = RemapMaliTopography(
+                test_case=self,
+                base_mesh_step=base_mesh_step,
+                name='remap_topo_unsmoothed',
+                mesh_name=mesh_name,
+                smoothing=False,
+                mali_ais_topo=mali_ais_topo,
                 ocean_includes_grounded=False)
 
-        self.add_step(remap_step)
+            self.add_step(unsmoothed_topo)
+
+            smoothed_topo = RemapMaliTopography(
+                test_case=self,
+                base_mesh_step=base_mesh_step,
+                name='remap_topo_unsmoothed',
+                mesh_name=mesh_name,
+                smoothing=False,
+                unsmoothed_topo=unsmoothed_topo,
+                mali_ais_topo=mali_ais_topo,
+                ocean_includes_grounded=False)
+
+            self.add_step(smoothed_topo)
 
         self.add_step(CullMeshStep(
             test_case=self, base_mesh_step=base_mesh_step,
             with_ice_shelf_cavities=self.with_ice_shelf_cavities,
-            remap_topography=remap_step))
+            unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
 
     def configure(self, config=None):
         """
@@ -174,14 +206,13 @@ def configure(self, config=None):
         if config is None:
             config = self.config
         config.add_from_package('compass.mesh', 'mesh.cfg', exception=True)
-        if 'remap_topography' in self.steps:
-            config.add_from_package('compass.ocean.mesh',
-                                    'remap_topography.cfg', exception=True)
+        config.add_from_package('compass.ocean.mesh',
+                                'remap_topography.cfg', exception=True)
 
-            if not self.high_res_topography:
-                config.add_from_package('compass.ocean.mesh',
-                                        'low_res_topography.cfg',
-                                        exception=True)
+        if not self.high_res_topography:
+            config.add_from_package('compass.ocean.mesh',
+                                    'low_res_topography.cfg',
+                                    exception=True)
 
         if self.mali_ais_topo is not None:
             package = 'compass.ocean.tests.global_ocean.mesh.' \
@@ -220,11 +251,7 @@ def configure(self, config=None):
                        'Antarctica')
 
         # a description of the bathymetry
-        if 'remap_topography' in self.steps:
-            description = config.get('remap_topography', 'description')
-        else:
-            description = 'Bathymetry is from GEBCO 2023, combined with ' \
-                          'BedMachine Antarctica v3 around Antarctica.'
+        description = config.get('remap_topography', 'description')
 
         config.set('global_ocean', 'bathy_description', description)