Merge pull request #141 from xylar/add-isomip-plus-base-and-culled-me…

…shes

Add steps for ISOMIP+ base and culled meshes

deploy/conda-dev-spec.template

@@ -35,7 +35,7 @@ otps={{ otps }}
 progressbar2
 pyamg>=4.2.2
 pyproj
-pyremap>=1.1.0,<2.0.0
+pyremap>=1.2.0,<2.0.0
 ruamel.yaml
 requests
 scipy>=1.8.0

polaris/ocean/__init__.py

@@ -5,6 +5,7 @@
 from polaris.ocean.tasks.inertial_gravity_wave import (
     add_inertial_gravity_wave_tasks,
 )
+from polaris.ocean.tasks.isomip_plus import add_isomip_plus_tasks
 from polaris.ocean.tasks.manufactured_solution import (
     add_manufactured_solution_tasks,
 )
@@ -26,6 +27,7 @@ def __init__(self):
         # planar: please keep these in alphabetical order
         add_baroclinic_channel_tasks(component=self)
         add_inertial_gravity_wave_tasks(component=self)
+        add_isomip_plus_tasks(component=self, mesh_type='planar')
         add_manufactured_solution_tasks(component=self)
 
         # single column
@@ -34,6 +36,7 @@ def __init__(self):
         # spherical: please keep these in alphabetical order
         add_cosine_bell_tasks(component=self)
         add_geostrophic_tasks(component=self)
+        add_isomip_plus_tasks(component=self, mesh_type='spherical')
         add_sphere_transport_tasks(component=self)
 
     def configure(self, config):

polaris/ocean/tasks/isomip_plus/__init__.py

@@ -0,0 +1,154 @@
+from typing import Dict, Union
+
+from polaris import Step
+from polaris.config import PolarisConfigParser
+from polaris.ocean.resolution import resolution_to_subdir
+from polaris.ocean.tasks.isomip_plus.isomip_plus_test import IsomipPlusTest
+from polaris.ocean.tasks.isomip_plus.mesh import (
+    CullMesh,
+    PlanarMesh,
+    SphericalMesh,
+)
+from polaris.ocean.tasks.isomip_plus.topo import TopoMap, TopoRemap, TopoScale
+
+
+def add_isomip_plus_tasks(component, mesh_type):
+    """
+    Add tasks for different baroclinic channel tests to the ocean component
+
+    component : polaris.ocean.Ocean
+        the ocean component that the tasks will be added to
+
+    mesh_type : {'planar', 'spherical'}
+        The type of mesh
+    """
+    planar = (mesh_type == 'planar')
+    for resolution in [4., 2., 1.]:
+        mesh_name = resolution_to_subdir(resolution)
+        resdir = f'{mesh_type}/isomip_plus/{mesh_name}'
+
+        filepath = f'{resdir}/isomip_plus_topo.cfg'
+        config = PolarisConfigParser(filepath=filepath)
+        if not planar:
+            config.add_from_package('polaris.mesh', 'spherical.cfg')
+            config.set('spherical_mesh', 'mpas_mesh_filename',
+                       'base_mesh_without_xy.nc')
+
+        config.add_from_package('polaris.ocean.tasks.isomip_plus',
+                                'isomip_plus.cfg')
+
+        config.add_from_package('polaris.ocean.tasks.isomip_plus',
+                                'isomip_plus_topo.cfg')
+
+        shared_steps = _get_shared_steps(
+            mesh_type, resolution, mesh_name, resdir, component, config)
+
+        for experiment in ['ocean0', 'ocean1', 'ocean2', 'ocean3', 'ocean4',
+                           'inception', 'wetting', 'drying']:
+            for vertical_coordinate in ['z-star']:
+                task = IsomipPlusTest(
+                    component=component,
+                    resdir=resdir,
+                    resolution=resolution,
+                    experiment=experiment,
+                    vertical_coordinate=vertical_coordinate,
+                    planar=planar,
+                    shared_steps=shared_steps[experiment])
+                component.add_task(task)
+
+
+def _get_shared_steps(mesh_type, resolution, mesh_name, resdir, component,
+                      config):
+    """ Get the shared steps for adding to tasks """
+
+    subdir = f'{resdir}/base_mesh'
+    base_mesh: Union[PlanarMesh, SphericalMesh, None] = None
+    if mesh_type == 'planar':
+        base_mesh = PlanarMesh(component=component,
+                               resolution=resolution,
+                               subdir=subdir,
+                               config=config)
+    else:
+        base_mesh = SphericalMesh(component=component,
+                                  cell_width=resolution,
+                                  subdir=subdir)
+        base_mesh.set_shared_config(config, link='isomip_plus_topo.cfg')
+
+    subdir = f'{resdir}/topo/map_base'
+    # we remap the topography onto the base mesh without smoothing and with
+    # the bilinear method because the conserve method doesn't work with
+    # periodic meshes
+    topo_map_base = TopoMap(component=component,
+                            name='topo_map_base',
+                            subdir=subdir,
+                            config=config,
+                            mesh_name=mesh_name,
+                            mesh_step=base_mesh,
+                            mesh_filename='base_mesh.nc',
+                            method='bilinear',
+                            smooth=False)
+
+    subdir = f'{resdir}/topo/remap_base'
+    topo_remap_base = TopoRemap(component=component,
+                                name='topo_remap_base',
+                                subdir=subdir,
+                                config=config,
+                                topo_map=topo_map_base,
+                                experiment='ocean1')
+
+    subdir = f'{resdir}/topo/cull_mesh'
+    cull_mesh = CullMesh(component=component,
+                         subdir=subdir,
+                         config=config,
+                         base_mesh=base_mesh,
+                         topo_remap=topo_remap_base)
+
+    subdir = f'{resdir}/topo/map_culled'
+    # we remap the topography onto the culled mesh with smoothing, which
+    # requires the conserve method
+    topo_map_culled = TopoMap(component=component,
+                              name='topo_map_culled',
+                              subdir=subdir,
+                              config=config,
+                              mesh_name=mesh_name,
+                              mesh_step=cull_mesh,
+                              mesh_filename='culled_mesh.nc',
+                              method='conserve',
+                              smooth=True)
+
+    topo_remap_culled: Dict[str, TopoRemap] = dict()
+    shared_steps: Dict[str, Dict[str, Step]] = dict()
+    for experiment in ['ocean1', 'ocean2', 'ocean3', 'ocean4']:
+        name = 'topo_remap_culled'
+        subdir = f'{resdir}/topo/remap_culled/{experiment}'
+        topo_remap_culled[experiment] = TopoRemap(component=component,
+                                                  name=name,
+                                                  subdir=subdir,
+                                                  config=config,
+                                                  topo_map=topo_map_culled,
+                                                  experiment=experiment)
+
+        shared_steps[experiment] = {
+            'base_mesh': base_mesh,
+            'topo/map_base': topo_map_base,
+            'topo/remap_base': topo_remap_base,
+            'topo/cull_mesh': cull_mesh,
+            'topo/map_culled': topo_map_culled,
+            'topo/remap_culled': topo_remap_culled[experiment],
+            'topo_final': topo_remap_culled[experiment]}
+
+    # ocean0 and ocean1 use the same topography
+    shared_steps['ocean0'] = shared_steps['ocean1']
+
+    for experiment in ['inception', 'wetting', 'drying']:
+        shared_steps[experiment] = dict(shared_steps['ocean1'])
+        subdir = f'{resdir}/topo/scale/{experiment}'
+        topo_scale = TopoScale(component=component,
+                               subdir=subdir,
+                               config=config,
+                               topo_remap=topo_remap_culled['ocean1'],
+                               experiment=experiment)
+        shared_steps[experiment]['topo/scale'] = topo_scale
+        shared_steps[experiment]['topo_final'] = topo_scale
+
+    return shared_steps

polaris/ocean/tasks/isomip_plus/isomip_plus.cfg

@@ -0,0 +1,5 @@
+# config options for ISOMIP+ test cases
+[isomip_plus]
+
+# the density of ice prescribed in ISOMIP+
+ice_density = 918

polaris/ocean/tasks/isomip_plus/isomip_plus_test.py

@@ -0,0 +1,92 @@
+from polaris import Task
+
+
+class IsomipPlusTest(Task):
+    """
+    An ISOMIP+ test case
+
+    Attributes
+    ----------
+    resolution : float
+        The horizontal resolution (km) of the test case
+
+    experiment : str
+        The ISOMIP+ experiment
+
+    vertical_coordinate : str
+            The type of vertical coordinate (``z-star``, ``z-level``, etc.)
+
+    tidal_forcing: bool
+        Whether the case has tidal forcing
+
+    thin_film: bool
+        Whether a thin film is present under land ice
+
+    planar : bool, optional
+        Whether the test case runs on a planar or a spherical mesh
+    """
+
+    def __init__(self, component, resdir, resolution, experiment,
+                 vertical_coordinate, planar, shared_steps,
+                 thin_film=False, tidal_forcing=False):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        component : polaris.Component
+            The component the task belongs to
+
+        resdir : str
+            The subdirectory in the component for ISOMIP+ experiments of the
+            given resolution
+
+        resolution : float
+            The horizontal resolution (km) of the test case
+
+        experiment : str
+            The ISOMIP+ experiment
+
+        vertical_coordinate : str
+            The type of vertical coordinate (``z-star``, ``z-level``, etc.)
+
+        planar : bool
+            Whether the test case runs on a planar or a spherical mesh
+
+        shared_steps : dict
+            The shared step for creating a topography mapping file from
+            the ISOMIP+ input data to the base mesh
+
+        thin_film: bool, optional
+            Whether the run includes a thin film below grounded ice
+
+        tidal_forcing: bool, optional
+            Whether the run includes a single-period tidal forcing
+        """  # noqa: E501
+        name = experiment
+        if tidal_forcing:
+            name = f'tidal_forcing_{name}'
+        if thin_film:
+            name = f'thin_film_{name}'
+
+        self.resolution = resolution
+        self.experiment = experiment
+        self.vertical_coordinate = vertical_coordinate
+        self.thin_film = thin_film
+        self.tidal_forcing = tidal_forcing
+        self.planar = planar
+        subdir = f'{resdir}/{vertical_coordinate}/{name}'
+        super().__init__(component=component, name=name, subdir=subdir)
+
+        for symlink, step in shared_steps.items():
+            if symlink == 'topo_final':
+                continue
+            self.add_step(step, symlink=symlink)
+
+    def configure(self):
+        """
+        Modify the configuration options for this test case.
+        """
+        config = self.config
+        config.add_from_package('polaris.ocean.tasks.isomip_plus',
+                                'isomip_plus.cfg')

polaris/ocean/tasks/isomip_plus/isomip_plus_topo.cfg

@@ -0,0 +1,43 @@
+# config options for ISOMIP+ meshes
+[isomip_plus_mesh]
+
+# latitude in degrees for origin of the mesh
+lat0 = -75.
+
+# size of the domain in km
+lx = 800
+ly = 80
+
+# a buffer in km around the domain that will be culled based on the topography
+buffer = 80
+
+# config options for ISOMIP+ topography
+[isomip_plus_topo]
+
+# minimum ocean fraction (i.e. fraction of bathymetry below sea level) of an
+# MPAS cell below which it will be culled from the mesh
+min_ocean_fraction = 0.5
+
+# the expansion factor used to smooth the topography
+expand_factor = 2.0
+
+# config options for ISOMIP+ topography scaling
+[isomip_plus_scaling]
+
+# simple thickening and thinning experiments that involve scaling the Ocean1
+# landIcePressure and landIceDraft over time
+#
+# "inception" reference dates
+inception_dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00
+# scaling at each date
+inception_scales = 0.0, 1.0, 1.0
+
+# "drying" reference dates
+drying_dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00
+# scaling at each date
+drying_scales = 1.0, 2.0, 2.0
+
+# "wetting" reference dates
+wetting_dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00
+# scaling at each date
+wetting_scales = 1.0, 0.0, 0.0

polaris/ocean/tasks/isomip_plus/mesh/__init__.py

@@ -0,0 +1,3 @@
+from polaris.ocean.tasks.isomip_plus.mesh.cull import CullMesh
+from polaris.ocean.tasks.isomip_plus.mesh.planar import PlanarMesh
+from polaris.ocean.tasks.isomip_plus.mesh.spherical import SphericalMesh