diff --git a/.github/workflows/build_workflow.yml b/.github/workflows/build_workflow.yml index ae88a9f7eb..ecd3cd4be7 100644 --- a/.github/workflows/build_workflow.yml +++ b/.github/workflows/build_workflow.yml @@ -36,7 +36,7 @@ jobs: python-version: "3.10" - if: ${{ steps.skip_check.outputs.should_skip != 'true' }} id: file_changes - uses: trilom/file-changes-action@v1.2.3 + uses: trilom/file-changes-action@v1.2.4 with: output: ' ' - if: ${{ steps.skip_check.outputs.should_skip != 'true' }} diff --git a/compass/landice/mesh.py b/compass/landice/mesh.py index a839525059..942405e346 100644 --- a/compass/landice/mesh.py +++ b/compass/landice/mesh.py @@ -1,7 +1,15 @@ import time import jigsawpy +import mpas_tools.io import numpy as np +import xarray +from geometric_features import FeatureCollection, GeometricFeatures +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from mpas_tools.mesh.conversion import convert, cull +from mpas_tools.mesh.creation import build_planar_mesh +from netCDF4 import Dataset def gridded_flood_fill(field, iStart=None, jStart=None): @@ -17,9 +25,11 @@ def gridded_flood_fill(field, iStart=None, jStart=None): field : numpy.ndarray Array from gridded dataset to use for flood-fill. Usually ice thickness. + iStart : int x index from which to start flood fill for field. Defaults to the center x coordinate. + jStart : int y index from which to start flood fill. Defaults to the center y coordinate. @@ -27,7 +37,7 @@ def gridded_flood_fill(field, iStart=None, jStart=None): Returns ------- flood_mask : numpy.ndarray - _mask calculated by the flood fill routine, + mask calculated by the flood fill routine, where cells connected to the ice sheet (or main feature) are 1 and everything else is 0. """ @@ -82,19 +92,23 @@ def set_rectangular_geom_points_and_edges(xmin, xmax, ymin, ymax): ---------- xmin : int or float Left-most x-coordinate in region to mesh + xmax : int or float Right-most x-coordinate in region to mesh + ymin : int or float Bottom-most y-coordinate in region to mesh + ymax : int or float Top-most y-coordinate in region to mesh Returns ------- geom_points : jigsawpy.jigsaw_msh_t.VERT2_t - xy node coordinates to pass to build_planar_mesh() + xy node coordinates to pass to ``build_planar_mesh()`` + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t - xy edge coordinates between nodes to pass to build_planar_mesh() + xy edge coordinates between nodes to pass to ``build_planar_mesh()`` """ geom_points = np.array([ # list of xy "node" coordinates @@ -114,45 +128,57 @@ def set_rectangular_geom_points_and_edges(xmin, xmax, ymin, ymax): return geom_points, geom_edges -def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, +def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None, dist_to_edge=None, dist_to_grounding_line=None, flood_fill_iStart=None, flood_fill_jStart=None): """ Set cell widths based on settings in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. - Requires the following options to be set in the given config section: - ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, - ``high_dist``, ``low_dist``,``cull_distance``, ``use_speed``, - ``use_dist_to_edge``, and ``use_dist_to_grounding_line``. Parameters ---------- - section : str - Section of the config file from which to read parameters + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + thk : numpy.ndarray Ice thickness field from gridded dataset, usually after trimming to flood fill mask + bed : numpy.ndarray Bed topography from gridded dataset + vx : numpy.ndarray, optional x-component of ice velocity from gridded dataset, usually after trimming to flood fill mask. Can be set to ``None`` if ``use_speed == False`` in config file. + vy : numpy.ndarray, optional y-component of ice velocity from gridded dataset, usually after trimming to flood fill mask. Can be set to ``None`` if ``use_speed == False`` in config file. + dist_to_edge : numpy.ndarray, optional Distance from each cell to ice edge, calculated in separate function. Can be set to ``None`` if ``use_dist_to_edge == False`` in config file and you do not want to set large ``cell_width`` where cells will be culled anyway, but this is not recommended. + dist_to_grounding_line : numpy.ndarray, optional Distance from each cell to grounding line, calculated in separate function. Can be set to ``None`` if ``use_dist_to_grounding_line == False`` in config file. + flood_fill_iStart : int, optional x-index location to start flood-fill when using bed topography + flood_fill_jStart : int, optional y-index location to start flood-fill when using bed topography @@ -164,19 +190,19 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, """ logger = self.logger - section = self.config[section] + section = self.config[section_name] # Get config inputs for cell spacing functions - min_spac = float(section.get('min_spac')) - max_spac = float(section.get('max_spac')) - high_log_speed = float(section.get('high_log_speed')) - low_log_speed = float(section.get('low_log_speed')) - high_dist = float(section.get('high_dist')) - low_dist = float(section.get('low_dist')) - high_dist_bed = float(section.get('high_dist_bed')) - low_dist_bed = float(section.get('low_dist_bed')) - low_bed = float(section.get('low_bed')) - high_bed = float(section.get('high_bed')) + min_spac = section.getfloat('min_spac') + max_spac = section.getfloat('max_spac') + high_log_speed = section.getfloat('high_log_speed') + low_log_speed = section.getfloat('low_log_speed') + high_dist = section.getfloat('high_dist') + low_dist = section.getfloat('low_dist') + high_dist_bed = section.getfloat('high_dist_bed') + low_dist_bed = section.getfloat('low_dist_bed') + low_bed = section.getfloat('low_bed') + high_bed = section.getfloat('high_bed') # convert km to m cull_distance = float(section.get('cull_distance')) * 1.e3 @@ -291,44 +317,59 @@ def set_cell_width(self, section, thk, bed=None, vx=None, vy=None, return cell_width -def get_dist_to_edge_and_GL(self, thk, topg, x, y, section, window_size=None): +def get_dist_to_edge_and_gl(self, thk, topg, x, y, + section_name, window_size=None): """ Calculate distance from each point to ice edge and grounding line, to be used in mesh density functions in :py:func:`compass.landice.mesh.set_cell_width()`. In future development, - this should be updated to use a faster package such as `scikit-fmm`. + this should be updated to use a faster package such as ``scikit-fmm``. Parameters ---------- thk : numpy.ndarray Ice thickness field from gridded dataset, usually after trimming to flood fill mask + topg : numpy.ndarray Bed topography field from gridded dataset + x : numpy.ndarray x coordinates from gridded dataset + y : numpy.ndarray y coordinates from gridded dataset - section : str - section of config file used to define mesh parameters + + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + window_size : int or float Size (in meters) of a search 'box' (one-directional) to use to calculate the distance from each cell to the ice margin. Bigger number makes search slower, but if too small, the transition zone could get truncated. We usually want this calculated as the - maximum of high_dist and high_dist_bed, but there may be cases in - which it is useful to set it manually. However, it should never be - smaller than either high_dist or high_dist_bed. + maximum of ``high_dist`` and ``high_dist_bed``, but there may be cases + in which it is useful to set it manually. However, it should never be + smaller than either ``high_dist`` or ``high_dist_bed``. Returns ------- dist_to_edge : numpy.ndarray Distance from each cell to the ice edge + dist_to_grounding_line : numpy.ndarray Distance from each cell to the grounding line """ logger = self.logger - section = self.config[section] + section = self.config[section_name] tic = time.time() high_dist = float(section.get('high_dist')) @@ -407,7 +448,299 @@ def get_dist_to_edge_and_GL(self, thk, topg, x, y, section, window_size=None): dist_to_grounding_line[i, j] = dist_to_here_grounding_line.min() toc = time.time() - logger.info('compass.landice.mesh.get_dist_to_edge_and_GL() took {:0.2f} ' + logger.info('compass.landice.mesh.get_dist_to_edge_and_gl() took {:0.2f} ' 'seconds'.format(toc - tic)) return dist_to_edge, dist_to_grounding_line + + +def build_cell_width(self, section_name, gridded_dataset, + flood_fill_start=[None, None]): + """ + Determine MPAS mesh cell size based on user-defined density function. + + Parameters + ---------- + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + + gridded_dataset : str + name of NetCDF file used to define cell spacing + + flood_fill_start : list of ints + ``i`` and ``j`` indices used to define starting location for flood + fill. Most cases will use ``[None, None]``, which will just start the + flood fill in the center of the gridded dataset. + + Returns + ------- + cell_width : numpy.ndarray + Desired width of MPAS cells based on mesh desnity functions to pass to + :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + + x1 : float + x coordinates from gridded dataset + + y1 : float + y coordinates from gridded dataset + + geom_points : jigsawpy.jigsaw_msh_t.VERT2_t + xy node coordinates to pass to ``build_planar_mesh()`` + + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t + xy edge coordinates between nodes to pass to ``build_planar_mesh()`` + + flood_mask : numpy.ndarray + mask calculated by the flood fill routine, + where cells connected to the ice sheet (or main feature) + are 1 and everything else is 0. + """ + + section = self.config[section_name] + # get needed fields from gridded dataset + f = Dataset(gridded_dataset, 'r') + f.set_auto_mask(False) # disable masked arrays + + x1 = f.variables['x1'][:] + y1 = f.variables['y1'][:] + thk = f.variables['thk'][0, :, :] + topg = f.variables['topg'][0, :, :] + vx = f.variables['vx'][0, :, :] + vy = f.variables['vy'][0, :, :] + + f.close() + + # Get bounds defined by user, or use bound of gridded dataset + bnds = [np.min(x1), np.max(x1), np.min(y1), np.max(y1)] + bnds_options = ['x_min', 'x_max', 'y_min', 'y_max'] + for index, option in enumerate(bnds_options): + bnd = section.get(option) + if bnd != 'None': + bnds[index] = float(bnd) + + geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds) + + # Remove ice not connected to the ice sheet. + flood_mask = gridded_flood_fill(thk) == 0 + thk[flood_mask] = 0.0 + vx[flood_mask] = 0.0 + vy[flood_mask] = 0.0 + + # Calculate distance from each grid point to ice edge + # and grounding line, for use in cell spacing functions. + distToEdge, distToGL = get_dist_to_edge_and_gl( + self, thk, topg, x1, + y1, section_name=section_name) + + # Set cell widths based on mesh parameters set in config file + cell_width = set_cell_width(self, section_name=section_name, + thk=thk, bed=topg, vx=vx, vy=vy, + dist_to_edge=distToEdge, + dist_to_grounding_line=distToGL, + flood_fill_iStart=flood_fill_start[0], + flood_fill_jStart=flood_fill_start[1]) + + return (cell_width.astype('float64'), x1.astype('float64'), + y1.astype('float64'), geom_points, geom_edges, flood_mask) + + +def build_mali_mesh(self, cell_width, x1, y1, geom_points, + geom_edges, mesh_name, section_name, + gridded_dataset, projection, geojson_file=None): + """ + Create the MALI mesh based on final cell widths determined by + :py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw and + MPAS-Tools functions. Culls the mesh based on config options, interpolates + all available fields from the gridded dataset to the MALI mesh using the + bilinear method, and marks domain boundaries as Dirichlet cells. + + Parameters + ---------- + cell_width : numpy.ndarray + Desired width of MPAS cells calculated by :py:func:`build_cell_width()` + based on mesh density functions define in :py:func:`set_cell_width()` + to pass to + :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + + x1 : float + x coordinates from gridded dataset + + y1 : float + y coordinates from gridded dataset + + geom_points : jigsawpy.jigsaw_msh_t.VERT2_t + xy node coordinates to pass to ``build_planar_mesh()`` + + geom_edges : jigsawpy.jigsaw_msh_t.EDGE2_t + xy edge coordinates between nodes to pass to ``build_planar_mesh()`` + + mesh_name : str + Filename to be used for final MALI NetCDF mesh file. + + section_name : str + Section of the config file from which to read parameters. The + following options to be set in the given config section: + ``levels``, ``x_min``, ``x_max``, ``y_min``, ``y_max``, + ``min_spac``, ``max_spac``, ``high_log_speed``, ``low_log_speed``, + ``high_dist``, ``low_dist``, ``high_dist_bed``, ``low_dist_bed``, + ``high_bed``, ``low_bed``, ``cull_distance``, ``use_speed``, + ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + See the Land-Ice Framework section of the Users or Developers guide + for more information about these options and their uses. + + gridded_dataset : str + Name of gridded dataset file to be used for interpolation to MALI mesh + + projection : str + Projection to be used for setting lat-long fields. + Likely ``'gis-gimp'`` or ``'ais-bedmap2'`` + + geojson_file : str + Name of geojson file that defines regional domain extent. + """ + + logger = self.logger + section = self.config[section_name] + + logger.info('calling build_planar_mesh') + build_planar_mesh(cell_width, x1, y1, geom_points, + geom_edges, logger=logger) + dsMesh = xarray.open_dataset('base_mesh.nc') + logger.info('culling mesh') + dsMesh = cull(dsMesh, logger=logger) + logger.info('converting to MPAS mesh') + dsMesh = convert(dsMesh, logger=logger) + logger.info('writing grid_converted.nc') + write_netcdf(dsMesh, 'grid_converted.nc') + levels = section.get('levels') + args = ['create_landice_grid_from_generic_MPAS_grid.py', + '-i', 'grid_converted.nc', + '-o', 'grid_preCull.nc', + '-l', levels, '-v', 'glimmer'] + + check_call(args, logger=logger) + + args = ['interpolate_to_mpasli_grid.py', '-s', + gridded_dataset, '-d', + 'grid_preCull.nc', '-m', 'b', '-t'] + + check_call(args, logger=logger) + + cullDistance = section.get('cull_distance') + if float(cullDistance) > 0.: + args = ['define_cullMask.py', '-f', + 'grid_preCull.nc', '-m', + 'distance', '-d', cullDistance] + + check_call(args, logger=logger) + else: + logger.info('cullDistance <= 0 in config file. ' + 'Will not cull by distance to margin. \n') + + if geojson_file is not None: + # This step is only necessary because the GeoJSON region + # is defined by lat-lon. + args = ['set_lat_lon_fields_in_planar_grid.py', '-f', + 'grid_preCull.nc', '-p', projection] + + check_call(args, logger=logger) + + args = ['MpasMaskCreator.x', 'grid_preCull.nc', + 'mask.nc', '-f', geojson_file] + + check_call(args, logger=logger) + + logger.info('culling to geojson file') + + dsMesh = xarray.open_dataset('grid_preCull.nc') + if geojson_file is not None: + mask = xarray.open_dataset('mask.nc') + else: + mask = None + + dsMesh = cull(dsMesh, dsInverse=mask, logger=logger) + write_netcdf(dsMesh, 'culled.nc') + + logger.info('Marking horns for culling') + args = ['mark_horns_for_culling.py', '-f', 'culled.nc'] + + check_call(args, logger=logger) + + logger.info('culling and converting') + dsMesh = xarray.open_dataset('culled.nc') + dsMesh = cull(dsMesh, logger=logger) + dsMesh = convert(dsMesh, logger=logger) + write_netcdf(dsMesh, 'dehorned.nc') + + args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', + 'dehorned.nc', '-o', + mesh_name, '-l', levels, '-v', 'glimmer', + '--beta', '--thermal', '--obs', '--diri'] + + check_call(args, logger=logger) + + args = ['interpolate_to_mpasli_grid.py', '-s', + gridded_dataset, '-d', mesh_name, '-m', 'b'] + + check_call(args, logger=logger) + + logger.info('Marking domain boundaries dirichlet') + args = ['mark_domain_boundaries_dirichlet.py', + '-f', mesh_name] + check_call(args, logger=logger) + + args = ['set_lat_lon_fields_in_planar_grid.py', '-f', + mesh_name, '-p', projection] + check_call(args, logger=logger) + + +def make_region_masks(self, mesh_filename, mask_filename, cores, tags): + """ + Create masks for ice-sheet subregions based on data + in ``MPAS-Dev/geometric_fatures``. + + Parameters + ---------- + mesh_filename : str + name of NetCDF mesh file for which to create region masks + + mask_filename : str + name of NetCDF file to contain region masks + + cores : int + number of processors used to create region masks + + tags : list of str + Groups of regions for which masks are to be defined + """ + + logger = self.logger + logger.info('creating region masks') + gf = GeometricFeatures() + fcMask = FeatureCollection() + + for tag in tags: + fc = gf.read(componentName='landice', objectType='region', + tags=[tag]) + fcMask.merge(fc) + + geojson_filename = 'regionMask.geojson' + fcMask.to_geojson(geojson_filename) + + args = ['compute_mpas_region_masks', + '-m', mesh_filename, + '-g', geojson_filename, + '-o', mask_filename, + '-t', 'cell', + '--process_count', f'{cores}', + '--format', mpas_tools.io.default_format, + '--engine', mpas_tools.io.default_engine] + check_call(args, logger=logger) diff --git a/compass/landice/tests/antarctica/mesh.py b/compass/landice/tests/antarctica/mesh.py index 5639d3c559..f36afe0979 100644 --- a/compass/landice/tests/antarctica/mesh.py +++ b/compass/landice/tests/antarctica/mesh.py @@ -1,20 +1,11 @@ from shutil import copyfile -import matplotlib.pyplot as plt -import mpas_tools import netCDF4 -import xarray -from geometric_features import FeatureCollection, GeometricFeatures -from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.mesh.creation import build_planar_mesh from compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, + build_cell_width, + build_mali_mesh, + make_region_masks, ) from compass.model import make_graph_file from compass.step import Step @@ -62,29 +53,13 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['antarctica'] + section_name = 'mesh' logger.info('calling build_cell_width') cell_width, x1, y1, geom_points, geom_edges, floodFillMask = \ - self.build_cell_width() - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'ais_8km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) + build_cell_width( + self, section_name=section_name, + gridded_dataset='antarctica_8km_2020_10_20.nc') # Apply floodFillMask to thickness field to help with culling copyfile('antarctica_8km_2020_10_20.nc', @@ -96,159 +71,25 @@ def run(self): gg.variables['vy'][0, :, :] *= floodFillMask gg.close() - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_8km_2020_10_20_floodFillMask.nc', '-d', - 'ais_8km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # Cull a certain distance from the ice margin - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'ais_8km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - dsMesh = xarray.open_dataset('ais_8km_preCull.nc') - dsMesh = cull(dsMesh, logger=logger) - write_netcdf(dsMesh, 'antarctica_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'antarctica_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('antarctica_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'antarctica_dehorned.nc') - - mesh_filename = self.mesh_filename - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'antarctica_dehorned.nc', '-o', - mesh_filename, '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_8km_2020_10_20.nc', - '-d', mesh_filename, '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', mesh_filename] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - mesh_filename, '-p', 'ais-bedmap2'] - check_call(args, logger=logger) + build_mali_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=self.mesh_filename, section_name=section_name, + gridded_dataset='antarctica_8km_2020_10_20_floodFillMask.nc', + projection='ais-bedmap2', geojson_file=None) logger.info('creating graph.info') - make_graph_file(mesh_filename=mesh_filename, + make_graph_file(mesh_filename=self.mesh_filename, graph_filename='graph.info') # create a region mask - mask_filename = f'{mesh_filename[:-3]}_imbie_regionMasks.nc' - self._make_region_masks(mesh_filename, mask_filename, - self.cpus_per_task, - tags=['EastAntarcticaIMBIE', - 'WestAntarcticaIMBIE', - 'AntarcticPeninsulaIMBIE']) - - mask_filename = f'{mesh_filename[:-3]}_ismip6_regionMasks.nc' - self._make_region_masks(mesh_filename, mask_filename, - self.cpus_per_task, - tags=['ISMIP6_Basin']) - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from AIS dataset - f = netCDF4.Dataset('antarctica_8km_2020_10_20.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Antarctica mesh. - xx0 = -3333500 - xx1 = 3330500 - yy0 = -3333500 - yy1 = 3330500 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL(self, thk, topg, x1, y1, - section='antarctica', - window_size=1.e5) - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section='antarctica', thk=thk, - vx=vx, vy=vy, dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges, floodMask) - - def _make_region_masks(self, mesh_filename, mask_filename, cores, tags): - logger = self.logger - gf = GeometricFeatures() - fcMask = FeatureCollection() - - for tag in tags: - fc = gf.read(componentName='landice', objectType='region', - tags=[tag]) - fc.plot('southpole') - plt.savefig(f'plot_basins_{tag}.png') - fcMask.merge(fc) - - geojson_filename = 'regionMask.geojson' - fcMask.to_geojson(geojson_filename) - - # these defaults may have been updated from config options -- pass them - # along to the subprocess - netcdf_format = mpas_tools.io.default_format - netcdf_engine = mpas_tools.io.default_engine - - args = ['compute_mpas_region_masks', - '-m', mesh_filename, - '-g', geojson_filename, - '-o', mask_filename, - '-t', 'cell', - '--process_count', f'{cores}', - '--format', netcdf_format, - '--engine', netcdf_engine] - check_call(args, logger=logger) + mask_filename = f'{self.mesh_filename[:-3]}_imbie_regionMasks.nc' + make_region_masks(self, self.mesh_filename, mask_filename, + self.cpus_per_task, + tags=['EastAntarcticaIMBIE', + 'WestAntarcticaIMBIE', + 'AntarcticPeninsulaIMBIE']) + + mask_filename = f'{self.mesh_filename[:-3]}_ismip6_regionMasks.nc' + make_region_masks(self, self.mesh_filename, mask_filename, + self.cpus_per_task, + tags=['ISMIP6_Basin']) diff --git a/compass/landice/tests/antarctica/antarctica.cfg b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg similarity index 71% rename from compass/landice/tests/antarctica/antarctica.cfg rename to compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg index 6bb0910815..bc9d62530d 100644 --- a/compass/landice/tests/antarctica/antarctica.cfg +++ b/compass/landice/tests/antarctica/mesh_gen/mesh_gen.cfg @@ -1,9 +1,15 @@ # config options for antarctica test cases -[antarctica] +[mesh] # number of levels in the mesh levels = 5 +# Bounds of Antarctic mesh +x_min = -3333500. +x_max = 3330500. +y_min = -3333500. +y_max = 3330500. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. @@ -22,10 +28,18 @@ low_log_speed = 0.75 high_dist = 1.e5 # distance to ice edge or grounding line within which cell spacing = min_spac (meters) low_dist = 1.e4 + +# These will not be applied unless use_bed = True. +# They are taken from the humboldt.mesh_gen test case +# and have not been evaluated for Antarctica. # distance at which bed topography has no effect high_dist_bed = 1.e5 # distance within which bed topography has maximum effect low_dist_bed = 7.5e4 +# Bed elev beneath which cell spacing is minimized +low_bed = 50.0 +# Bed elev above which cell spacing is maximized +high_bed = 100.0 # mesh density functions use_speed = True diff --git a/compass/landice/tests/greenland/__init__.py b/compass/landice/tests/greenland/__init__.py index 881438431d..500338d5ff 100644 --- a/compass/landice/tests/greenland/__init__.py +++ b/compass/landice/tests/greenland/__init__.py @@ -1,8 +1,10 @@ -from compass.testgroup import TestGroup -from compass.landice.tests.greenland.smoke_test import SmokeTest -from compass.landice.tests.greenland.decomposition_test import DecompositionTest +from compass.landice.tests.greenland.decomposition_test import ( + DecompositionTest, +) +from compass.landice.tests.greenland.mesh_gen import MeshGen from compass.landice.tests.greenland.restart_test import RestartTest -from compass.landice.tests.greenland.high_res_mesh import HighResMesh +from compass.landice.tests.greenland.smoke_test import SmokeTest +from compass.testgroup import TestGroup class Greenland(TestGroup): @@ -27,4 +29,4 @@ def __init__(self, mpas_core): RestartTest(test_group=self, velo_solver=velo_solver)) self.add_test_case( - HighResMesh(test_group=self)) + MeshGen(test_group=self)) diff --git a/compass/landice/tests/greenland/mesh.py b/compass/landice/tests/greenland/mesh.py index d3851b57c2..bb7dc702e9 100644 --- a/compass/landice/tests/greenland/mesh.py +++ b/compass/landice/tests/greenland/mesh.py @@ -1,16 +1,7 @@ -import netCDF4 -import numpy as np -import xarray -from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.mesh.creation import build_planar_mesh - from compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, + build_cell_width, + build_mali_mesh, + make_region_masks, ) from compass.model import make_graph_file from compass.step import Step @@ -56,141 +47,35 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_GIS_mesh'] - - logger.info('calling build_cell_wdith') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - cullDistance = section.get('cull_distance') - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - dsMesh = cull(dsMesh, logger=logger) - write_netcdf(dsMesh, 'greenland_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'greenland_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('greenland_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'greenland_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'greenland_dehorned.nc', '-o', - 'GIS.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'GIS.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'GIS.nc'] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'GIS.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + mesh_name = 'GIS.nc' + section_name = 'mesh' + + logger.info('calling build_cell_width') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name=section_name, + gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc', + flood_fill_start=[100, 700]) + + build_mali_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa + projection='gis-gimp', geojson_file=None) logger.info('creating graph.info') - make_graph_file(mesh_filename='GIS.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes the definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_2km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the GIS mesh. - xx0 = np.min(x1) - xx1 = np.max(x1) - yy0 = np.min(y1) - yy1 = np.max(y1) - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( - self, thk, topg, x1, y1, section='high_res_GIS_mesh') - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file. - # Start flood-fill for bed topography in the ocean at a point - # with bed < low_bed. This is often necessary to remove pockets - # of high resolution that are no longer connected to the rest of - # the high resolution part of the mesh. These indices are specific - # to the mesh and the gridded data product used to set cell spacing. - cell_width = set_cell_width(self, section='high_res_GIS_mesh', thk=thk, - bed=topg, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL, - flood_fill_iStart=100, - flood_fill_jStart=700) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) + # create region masks + mask_filename = f'{mesh_name[:-3]}_regionMasks.nc' + make_region_masks(self, mesh_name, mask_filename, + self.cpus_per_task, + tags=['eastCentralGreenland', + 'northEastGreenland', + 'northGreenland', + 'northWestGreenland', + 'southEastGreenland', + 'southGreenland', + 'southWestGreenland', + 'westCentralGreenland']) diff --git a/compass/landice/tests/greenland/high_res_mesh/__init__.py b/compass/landice/tests/greenland/mesh_gen/__init__.py similarity index 92% rename from compass/landice/tests/greenland/high_res_mesh/__init__.py rename to compass/landice/tests/greenland/mesh_gen/__init__.py index 8b09f9a13d..b8608994e6 100644 --- a/compass/landice/tests/greenland/high_res_mesh/__init__.py +++ b/compass/landice/tests/greenland/mesh_gen/__init__.py @@ -1,8 +1,8 @@ -from compass.testcase import TestCase from compass.landice.tests.greenland.mesh import Mesh +from compass.testcase import TestCase -class HighResMesh(TestCase): +class MeshGen(TestCase): """ The high resolution test case for the greenland test group simply creates the mesh and initial condition. @@ -18,7 +18,7 @@ def __init__(self, test_group): test_group : compass.landice.tests.greenland.Greenland The test group that this test case belongs to """ - name = 'high_res_mesh' + name = 'mesh_gen' subdir = name super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg b/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg similarity index 85% rename from compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg rename to compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg index 8ad755ba30..cea3358b46 100644 --- a/compass/landice/tests/greenland/high_res_mesh/high_res_mesh.cfg +++ b/compass/landice/tests/greenland/mesh_gen/mesh_gen.cfg @@ -1,9 +1,17 @@ # config options for high_res_mesh test case -[high_res_GIS_mesh] +[mesh] # number of levels in the mesh levels = 10 +# Bounds of GIS mesh. If any bound is set +# to None, the mesh will use the full extent +# of the gridded dataset. +x_min = None +x_max = None +y_min = None +y_max = None + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/humboldt/mesh.py b/compass/landice/tests/humboldt/mesh.py index 5d725aaa54..8c168f8c80 100644 --- a/compass/landice/tests/humboldt/mesh.py +++ b/compass/landice/tests/humboldt/mesh.py @@ -1,16 +1,4 @@ -import netCDF4 -import xarray -from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.mesh.creation import build_planar_mesh - -from compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, -) +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file from compass.step import Step @@ -39,7 +27,7 @@ def __init__(self, test_case): super().__init__(test_case=test_case, name='mesh') self.add_output_file(filename='graph.info') - self.add_output_file(filename='Humboldt_1to10km.nc') + self.add_output_file(filename='Humboldt.nc') self.add_input_file( filename='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', target='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', @@ -59,163 +47,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['humboldt_mesh'] - - logger.info('calling build_cell_wdith') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - # If no number of levels specified in config file, use 10 - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - # This step uses a subset of the whole Greenland dataset trimmed to - # the region around Humboldt Glacier, to speed up interpolation. - # This could also be replaced with the full Greenland Ice Sheet - # dataset. - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'gis_1km_preCull.nc', '-p', 'gis-gimp'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'gis_1km_preCull.nc', - 'humboldt_mask.nc', '-f', 'Humboldt.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - humboldtMask = xarray.open_dataset('humboldt_mask.nc') - dsMesh = cull(dsMesh, dsInverse=humboldtMask, logger=logger) - write_netcdf(dsMesh, 'humboldt_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'humboldt_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('humboldt_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'humboldt_dehorned.nc') + section_name = 'mesh' + mesh_name = 'Humboldt.nc' - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'humboldt_dehorned.nc', '-o', - 'Humboldt_1to10km.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] + logger.info('calling build_cell_width') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name=section_name, + gridded_dataset='greenland_2km_2020_04_20.epsg3413.nc') - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'Humboldt_1to10km.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Humboldt_1to10km.nc'] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'Humboldt_1to10km.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + build_mali_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='humboldt_1km_2020_04_20.epsg3413.icesheetonly.nc', + projection='gis-gimp', geojson_file='Humboldt.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Humboldt_1to10km.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. In the future, this function - and its components will likely be separated into separate generalized - functions to be reusable by multiple test groups. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_2km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Humboldt Glacier mesh. - xx0 = -630000 - xx1 = 84000 - yy0 = -1560000 - yy1 = -860000 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( - self, thk, topg, x1, - y1, section='humboldt_mesh') - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section='humboldt_mesh', - thk=thk, bed=topg, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/humboldt/humboldt.cfg b/compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg similarity index 90% rename from compass/landice/tests/humboldt/humboldt.cfg rename to compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg index 9bac801e01..6927514134 100644 --- a/compass/landice/tests/humboldt/humboldt.cfg +++ b/compass/landice/tests/humboldt/mesh_gen/mesh_gen.cfg @@ -1,9 +1,15 @@ # config options for humboldt test cases -[humboldt_mesh] +[mesh] # number of levels in the mesh levels = 10 +# Bounds of Humboldt domain +x_min = -630000. +x_max = 84000. +y_min = -1560000. +y_max = -860000. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/kangerlussuaq/mesh.py b/compass/landice/tests/kangerlussuaq/mesh.py index 8cd4b3c8a0..a142c34f79 100644 --- a/compass/landice/tests/kangerlussuaq/mesh.py +++ b/compass/landice/tests/kangerlussuaq/mesh.py @@ -1,16 +1,4 @@ -import netCDF4 -import xarray -from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.mesh.creation import build_planar_mesh - -from compass.landice.mesh import ( - get_dist_to_edge_and_GL, - gridded_flood_fill, - set_cell_width, - set_rectangular_geom_points_and_edges, -) +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file from compass.step import Step @@ -60,158 +48,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_Kangerlussuaq_mesh'] - - logger.info('calling build_cell_wdith') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'gis_1km_preCull.nc', '-p', 'gis-gimp'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'gis_1km_preCull.nc', - 'kangerlussuaq_mask.nc', '-f', 'Kangerlussuaq.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - kangerMask = xarray.open_dataset('kangerlussuaq_mask.nc') - dsMesh = cull(dsMesh, dsInverse=kangerMask, logger=logger) - write_netcdf(dsMesh, 'kangerlussuaq_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'kangerlussuaq_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('kangerlussuaq_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'kangerlussuaq_dehorned.nc') + mesh_name = 'Kangerlussuaq.nc' + section_name = 'mesh' - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'kangerlussuaq_dehorned.nc', '-o', - 'Kangerlussuaq.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] + logger.info('calling build_cell_width') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name=section_name, + gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'Kangerlussuaq.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Kangerlussuaq.nc'] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'Kangerlussuaq.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + build_mali_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa + projection='gis-gimp', geojson_file='Kangerlussuaq.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Kangerlussuaq.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_8km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Kangerlussuaq Glacier mesh. - xx0 = 215000 - xx1 = 544000 - yy0 = -2370000 - yy1 = -2070000 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL( - self, thk, topg, x1, - y1, section='high_res_Kangerlussuaq_mesh') - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, - section='high_res_Kangerlussuaq_mesh', - thk=thk, bed=topg, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=distToGL) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg b/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg index 7db6fd63a6..72c5a9da2f 100644 --- a/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/kangerlussuaq/mesh_gen/mesh_gen.cfg @@ -1,9 +1,15 @@ # config options for kangerlussuaq test cases -[high_res_Kangerlussuaq_mesh] +[mesh] # number of levels in the mesh levels = 10 +# Bounds of Kangerlussuaq domain +x_min = 215000. +x_max = 544000. +y_min = -2370000. +y_max = -2070000. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. diff --git a/compass/landice/tests/koge_bugt_s/mesh.py b/compass/landice/tests/koge_bugt_s/mesh.py index 9a6fa5442b..4a41fdb36b 100644 --- a/compass/landice/tests/koge_bugt_s/mesh.py +++ b/compass/landice/tests/koge_bugt_s/mesh.py @@ -1,19 +1,6 @@ -import numpy as np -import netCDF4 -import xarray -from matplotlib import pyplot as plt - -from mpas_tools.mesh.creation import build_planar_mesh -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.planar_hex import make_planar_hex_mesh -from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call - -from compass.step import Step +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file -from compass.landice.mesh import gridded_flood_fill, \ - set_rectangular_geom_points_and_edges, \ - set_cell_width, get_dist_to_edge_and_GL +from compass.step import Step class Mesh(Step): @@ -43,9 +30,9 @@ def __init__(self, test_case): self.add_output_file(filename='graph.info') self.add_output_file(filename='Koge_Bugt_S.nc') self.add_input_file( - filename='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - target='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - database='') + filename='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', + target='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', + database='') self.add_input_file(filename='Koge_Bugt_S.geojson', package='compass.landice.tests.koge_bugt_s', target='Koge_Bugt_S.geojson', @@ -61,157 +48,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_KogeBugtS_mesh'] + mesh_name = 'Koge_Bugt_S.nc' + section_name = 'mesh' logger.info('calling build_cell_width') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - dsMesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - dsMesh = cull(dsMesh, logger=logger) - logger.info('converting to MPAS mesh') - dsMesh = convert(dsMesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(dsMesh, 'grid_converted.nc') - levels = section.get('levels') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'gis_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', '-d', - 'gis_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'gis_1km_preCull.nc', '-m', - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'gis_1km_preCull.nc', '-p', 'gis-gimp'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'gis_1km_preCull.nc', - 'koge_bugt_s_mask.nc', '-f', 'Koge_Bugt_S.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - dsMesh = xarray.open_dataset('gis_1km_preCull.nc') - kangerMask = xarray.open_dataset('koge_bugt_s_mask.nc') - dsMesh = cull(dsMesh, dsInverse=kangerMask, logger=logger) - write_netcdf(dsMesh, 'koge_bugt_s_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'koge_bugt_s_culled.nc'] - check_call(args, logger=logger) - - logger.info('culling and converting') - dsMesh = xarray.open_dataset('koge_bugt_s_culled.nc') - dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) - write_netcdf(dsMesh, 'koge_bugt_s_dehorned.nc') + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name=section_name, + gridded_dataset='greenland_8km_2020_04_20.epsg3413.nc') - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'koge_bugt_s_dehorned.nc', '-o', - 'Koge_Bugt_S.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', - '-d', 'Koge_Bugt_S.nc', '-m', 'b', '-t'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Koge_Bugt_S.nc'] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'Koge_Bugt_S.nc', '-p', 'gis-gimp'] - check_call(args, logger=logger) + build_mali_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='greenland_1km_2020_04_20.epsg3413.icesheetonly.nc', # noqa + projection='gis-gimp', geojson_file='Koge_Bugt_S.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Koge_Bugt_S.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. - """ - # get needed fields from GIS dataset - f = netCDF4.Dataset('greenland_8km_2020_04_20.epsg3413.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Koge_Bugt_S Glacier mesh. - xx0 = 0 - xx1 = 225000 - yy0 = -2830000 - yy1 = -2700000 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - floodMask = gridded_flood_fill(thk) - thk[floodMask == 0] = 0.0 - vx[floodMask == 0] = 0.0 - vy[floodMask == 0] = 0.0 - - # Calculate distance from each grid point to ice edge - # and grounding line, for use in cell spacing functions. - distToEdge, distToGL = get_dist_to_edge_and_GL(self, thk, topg, x1, - y1, window_size=1.e5) - # optional - plot distance calculation - # plt.pcolor(distToEdge/1000.0); plt.colorbar(); plt.show() - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, - section='high_res_KogeBugtS_mesh', - thk=thk, vx=vx, vy=vy, - dist_to_edge=distToEdge, - dist_to_grounding_line=None) - # plt.pcolor(cell_width); plt.colorbar(); plt.show() - - return (cell_width.astype('float64'), x1.astype('float64'), - y1.astype('float64'), geom_points, geom_edges) diff --git a/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg b/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg index 4631a34063..770432a6d8 100644 --- a/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg +++ b/compass/landice/tests/koge_bugt_s/mesh_gen/mesh_gen.cfg @@ -1,9 +1,15 @@ # config options for koge_bugt_s test cases -[high_res_KogeBugtS_mesh] +[mesh] # number of levels in the mesh levels = 10 +# Bounds of Koge Bugt S mesh +x_min = 0. +x_max = 225000. +y_min = -2830000. +y_max = -2700000. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. @@ -23,7 +29,18 @@ high_dist = 1.e5 # distance within which cell spacing = min_spac (meters) low_dist = 1.e4 +# mesh density parameters used if use_bed = True +# distance at which bed topography has no effect +high_dist_bed = 1.e5 +# distance within which bed topography has maximum effect +low_dist_bed = 5.e4 +# Bed elev beneath which cell spacing is minimized +low_bed = 50.0 +# Bed elev above which cell spacing is maximized +high_bed = 100.0 + # mesh density functions use_speed = True use_dist_to_grounding_line = False use_dist_to_edge = True +use_bed = False diff --git a/compass/landice/tests/thwaites/__init__.py b/compass/landice/tests/thwaites/__init__.py index 74343beda3..27a764eaa9 100644 --- a/compass/landice/tests/thwaites/__init__.py +++ b/compass/landice/tests/thwaites/__init__.py @@ -1,7 +1,7 @@ -from compass.testgroup import TestGroup from compass.landice.tests.thwaites.decomposition_test import DecompositionTest +from compass.landice.tests.thwaites.mesh_gen import MeshGen from compass.landice.tests.thwaites.restart_test import RestartTest -from compass.landice.tests.thwaites.high_res_mesh import HighResMesh +from compass.testgroup import TestGroup class Thwaites(TestGroup): @@ -20,4 +20,4 @@ def __init__(self, mpas_core): self.add_test_case(RestartTest(test_group=self)) - self.add_test_case(HighResMesh(test_group=self)) + self.add_test_case(MeshGen(test_group=self)) diff --git a/compass/landice/tests/thwaites/mesh.py b/compass/landice/tests/thwaites/mesh.py index e0cb2da985..863bbd6796 100644 --- a/compass/landice/tests/thwaites/mesh.py +++ b/compass/landice/tests/thwaites/mesh.py @@ -1,16 +1,6 @@ -import netCDF4 -import xarray - -from mpas_tools.mesh.creation import build_planar_mesh -from mpas_tools.mesh.conversion import convert, cull -from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call - -from compass.step import Step +from compass.landice.mesh import build_cell_width, build_mali_mesh from compass.model import make_graph_file -from compass.landice.mesh import gridded_flood_fill, \ - set_rectangular_geom_points_and_edges, \ - set_cell_width, get_dist_to_edge_and_GL +from compass.step import Step class Mesh(Step): @@ -47,161 +37,21 @@ def run(self): Run this step of the test case """ logger = self.logger - config = self.config - section = config['high_res_mesh'] + mesh_name = 'Thwaites.nc' + section_name = 'mesh' logger.info('calling build_cell_width') - cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width() - logger.info('calling build_planar_mesh') - build_planar_mesh(cell_width, x1, y1, geom_points, - geom_edges, logger=logger) - ds_mesh = xarray.open_dataset('base_mesh.nc') - logger.info('culling mesh') - ds_mesh = cull(ds_mesh, logger=logger) - logger.info('converting to MPAS mesh') - ds_mesh = convert(ds_mesh, logger=logger) - logger.info('writing grid_converted.nc') - write_netcdf(ds_mesh, 'grid_converted.nc') - - levels = section.get('levels') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', - '-i', 'grid_converted.nc', - '-o', 'ase_1km_preCull.nc', - '-l', levels, '-v', 'glimmer'] - check_call(args, logger=logger) - - # This step uses a subset of the whole Antarctica dataset trimmed to - # the Amundsen Sean Embayment, to speed up interpolation. - # This could also be replaced with the full Antarctic Ice Sheet - # dataset. - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_1km_2020_10_20_ASE.nc', '-d', - 'ase_1km_preCull.nc', '-m', 'b', '-t'] - - check_call(args, logger=logger) - - # This step is only necessary if you wish to cull a certain - # distance from the ice margin, within the bounds defined by - # the GeoJSON file. - cullDistance = section.get('cull_distance') - if float(cullDistance) > 0.: - logger.info('calling define_cullMask.py') - args = ['define_cullMask.py', '-f', - 'ase_1km_preCull.nc', '-m' - 'distance', '-d', cullDistance] - - check_call(args, logger=logger) - else: - logger.info('cullDistance <= 0 in config file. ' - 'Will not cull by distance to margin. \n') - - # This step is only necessary because the GeoJSON region - # is defined by lat-lon. - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'ase_1km_preCull.nc', '-p', 'ais-bedmap2'] - - check_call(args, logger=logger) - - logger.info('calling MpasMaskCreator.x') - args = ['MpasMaskCreator.x', 'ase_1km_preCull.nc', - 'thwaites_mask.nc', '-f', 'thwaites_minimal.geojson'] - - check_call(args, logger=logger) - - logger.info('culling to geojson file') - ds_mesh = xarray.open_dataset('ase_1km_preCull.nc') - thwaitesMask = xarray.open_dataset('thwaites_mask.nc') - ds_mesh = cull(ds_mesh, dsInverse=thwaitesMask, logger=logger) - write_netcdf(ds_mesh, 'thwaites_culled.nc') - - logger.info('Marking horns for culling') - args = ['mark_horns_for_culling.py', '-f', 'thwaites_culled.nc'] - check_call(args, logger=logger) + cell_width, x1, y1, geom_points, geom_edges, floodMask = \ + build_cell_width( + self, section_name=section_name, + gridded_dataset='antarctica_8km_2020_10_20.nc') - logger.info('culling and converting') - ds_mesh = xarray.open_dataset('thwaites_culled.nc') - ds_mesh = cull(ds_mesh, logger=logger) - ds_mesh = convert(ds_mesh, logger=logger) - write_netcdf(ds_mesh, 'thwaites_dehorned.nc') - - logger.info('calling create_landice_grid_from_generic_MPAS_grid.py') - args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', - 'thwaites_dehorned.nc', '-o', - 'Thwaites.nc', '-l', levels, '-v', 'glimmer', - '--beta', '--thermal', '--obs', '--diri'] - - check_call(args, logger=logger) - - logger.info('calling interpolate_to_mpasli_grid.py') - args = ['interpolate_to_mpasli_grid.py', '-s', - 'antarctica_1km_2020_10_20_ASE.nc', - '-d', 'Thwaites.nc', '-m', 'b'] - check_call(args, logger=logger) - - logger.info('Marking domain boundaries dirichlet') - args = ['mark_domain_boundaries_dirichlet.py', - '-f', 'Thwaites.nc'] - check_call(args, logger=logger) - - logger.info('calling set_lat_lon_fields_in_planar_grid.py') - args = ['set_lat_lon_fields_in_planar_grid.py', '-f', - 'Thwaites.nc', '-p', 'ais-bedmap2'] - check_call(args, logger=logger) + build_mali_mesh( + self, cell_width, x1, y1, geom_points, geom_edges, + mesh_name=mesh_name, section_name=section_name, + gridded_dataset='antarctica_1km_2020_10_20_ASE.nc', + projection='ais-bedmap2', geojson_file='thwaites_minimal.geojson') logger.info('creating graph.info') - make_graph_file(mesh_filename='Thwaites.nc', + make_graph_file(mesh_filename=mesh_name, graph_filename='graph.info') - - def build_cell_width(self): - """ - Determine MPAS mesh cell size based on user-defined density function. - - This includes hard-coded definition of the extent of the regional - mesh and user-defined mesh density functions based on observed flow - speed and distance to the ice margin. In the future, this function - and its components will likely be separated into separate generalized - functions to be reusable by multiple test groups. - """ - # get needed fields from Antarctica dataset - f = netCDF4.Dataset('antarctica_8km_2020_10_20.nc', 'r') - f.set_auto_mask(False) # disable masked arrays - config = self.config - - x1 = f.variables['x1'][:] - y1 = f.variables['y1'][:] - thk = f.variables['thk'][0, :, :] - topg = f.variables['topg'][0, :, :] - vx = f.variables['vx'][0, :, :] - vy = f.variables['vy'][0, :, :] - - # Define extent of region to mesh. - # These coords are specific to the Thwaites Glacier mesh. - xx0 = -1864434 - xx1 = -975432 - yy0 = -901349 - yy1 = 0 - geom_points, geom_edges = set_rectangular_geom_points_and_edges( - xx0, xx1, yy0, yy1) - - # Remove ice not connected to the ice sheet. - flood_mask = gridded_flood_fill(thk) - thk[flood_mask == 0] = 0.0 - vx[flood_mask == 0] = 0.0 - vy[flood_mask == 0] = 0.0 - - # Calculate distances to ice edge and grounding line - dist_to_edge, dist_to_GL = get_dist_to_edge_and_GL(self, thk, - topg, x1, y1, - window_size=1.e5) - - # Set cell widths based on mesh parameters set in config file - cell_width = set_cell_width(self, section='high_res_mesh', thk=thk, - vx=vx, vy=vy, dist_to_edge=dist_to_edge, - dist_to_grounding_line=dist_to_GL) - - return cell_width.astype('float64'), x1.astype('float64'), \ - y1.astype('float64'), geom_points, geom_edges diff --git a/compass/landice/tests/thwaites/high_res_mesh/__init__.py b/compass/landice/tests/thwaites/mesh_gen/__init__.py similarity index 92% rename from compass/landice/tests/thwaites/high_res_mesh/__init__.py rename to compass/landice/tests/thwaites/mesh_gen/__init__.py index 8a44ee7f62..184a64642b 100644 --- a/compass/landice/tests/thwaites/high_res_mesh/__init__.py +++ b/compass/landice/tests/thwaites/mesh_gen/__init__.py @@ -1,8 +1,8 @@ -from compass.testcase import TestCase from compass.landice.tests.thwaites.mesh import Mesh +from compass.testcase import TestCase -class HighResMesh(TestCase): +class MeshGen(TestCase): """ The high resolution test case for the thwaites test group simply creates the mesh and initial condition. @@ -18,7 +18,7 @@ def __init__(self, test_group): test_group : compass.landice.tests.thwaites.Thwaites The test group that this test case belongs to """ - name = 'high_res_mesh' + name = 'mesh_gen' subdir = name super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg b/compass/landice/tests/thwaites/mesh_gen/mesh_gen.cfg similarity index 56% rename from compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg rename to compass/landice/tests/thwaites/mesh_gen/mesh_gen.cfg index dcb2ac81a0..98ae96049e 100644 --- a/compass/landice/tests/thwaites/high_res_mesh/high_res_mesh.cfg +++ b/compass/landice/tests/thwaites/mesh_gen/mesh_gen.cfg @@ -1,9 +1,15 @@ # config options for high_res_mesh test case -[high_res_mesh] +[mesh] # number of levels in the mesh levels = 10 +# Bounds of Thwaites regional mesh +x_min = -1864434. +x_max = -975432. +y_min = -901349. +y_max = 0. + # distance from ice margin to cull (km). # Set to a value <= 0 if you do not want # to cull based on distance from margin. @@ -23,7 +29,20 @@ high_dist = 1.e5 # distance within which cell spacing = min_spac low_dist = 5.e4 +# mesh density parameters used if use_bed = True +# These settings are taken from the Humboldt mesh +# and have not yet been evaluated for Thwaites. +# distance at which bed topography has no effect +high_dist_bed = 1.e5 +# distance within which bed topography has maximum effect +low_dist_bed = 5.e4 +# Bed elev beneath which cell spacing is minimized +low_bed = 50.0 +# Bed elev above which cell spacing is maximized +high_bed = 100.0 + # mesh density functions use_speed = True use_dist_to_grounding_line = True use_dist_to_edge = False +use_bed = False diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index bddb48c64a..9a51322d0a 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -203,10 +203,9 @@ greenland mesh.Mesh mesh.Mesh.run - mesh.Mesh.build_cell_width - high_res_mesh.HighResMesh - high_res_mesh.HighResMesh.run + mesh_gen.MeshGen + mesh_gen.MeshGen.run humboldt ~~~~~~~~ @@ -400,6 +399,12 @@ thwaites run_model.RunModel.setup run_model.RunModel.run + mesh_gen.MeshGen + mesh_gen.MeshGen.run + + mesh.Mesh + mesh.Mesh.run + Landice Framework ^^^^^^^^^^^^^^^^^ @@ -413,4 +418,7 @@ Landice Framework mesh.gridded_flood_fill mesh.set_rectangular_geom_points_and_edges mesh.set_cell_width - mesh.get_dist_to_edge_and_GL + mesh.get_dist_to_edge_and_gl + mesh.build_cell_width + mesh.build_mali_mesh + mesh.make_region_masks diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index 2bb6433ca7..f8c0fc858c 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -10,3 +10,98 @@ extrapolate The landice framework module ``compass/landice/extrapolate.py`` provides a function for extrapolating variables into undefined regions. It is copied from a similar script in MPAS-Tools. + +mesh +~~~~ + +The landice framework module :py:mod:`compass.landice.mesh` provides functions +used by all ``mesh_gen`` tests cases, which currently exist within the +``antarctica``, ``greenland``, ``humboldt``, ``kangerlussuaq``, ``koge_bugt_s``, +and ``thwaites`` test groups. These functions include: + +:py:func:`compass.landice.mesh.gridded_flood_fill()` applies a flood-fill algorithm +to the gridded dataset in order to separate the ice sheet from peripheral ice. + +:py:func:`compass.landice.mesh.set_rectangular_geom_points_and_edges()` sets node +and edge coordinates to pass to py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. + +:py:func:`compass.landice.mesh.set_cell_width()` sets cell widths based on settings +in config file to pass to :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. +Requires the following options to be set in the given config section: ``min_spac``, +``max_spac``, ``high_log_speed``, ``low_log_speed``, ``high_dist``, ``low_dist``, +``high_dist_bed``, ``low_dist_bed``, ``high_bed``, ``low_bed``, ``cull_distance``, +``use_speed``, ``use_dist_to_edge``, ``use_dist_to_grounding_line``, and ``use_bed``. + +:py:func:`compass.landice.mesh.get_dist_to_edge_and_gl()` calculates distance from +each point to ice edge and grounding line, to be used in mesh density functions in +:py:func:`compass.landice.mesh.set_cell_width()`. In future development, +this should be updated to use a faster package such as `scikit-fmm`. + +:py:func:`compass.landice.mesh.build_cell_width()` determine final MPAS mesh cell sizes +using desired cell widths calculated by py:func:`compass.landice.mesh.set_cell_width()`, +based on user-defined density functions and config options. + +:py:func:`compass.landice.mesh.build_mali_mesh()` creates the MALI mesh based on final +cell widths determined by py:func:`compass.landice.mesh.build_cell_width()`, using Jigsaw +and MPAS-Tools functions. Culls the mesh based on config options, interpolates +all available fields from the gridded dataset to the MALI mesh using the bilinear +method, and marks domain boundaries as Dirichlet cells. + +:py:func:`compass.landice.mesh.make_region_masks()` creates region masks using regions +defined in Geometric Features repository. It is only used by the ``antarctica`` +and ``greenland`` test cases. + +The following config options should be defined for all ``mesh_gen`` test cases (although +not necessarily with the same values shown here, which are the defaults for the 1–10km +Humboldt mesh): + +.. code-block:: cfg + + # config options for humboldt test cases + [mesh] + + # number of levels in the mesh + levels = 10 + + # Bounds of Humboldt domain. If you want the extent + # of the gridded dataset to determine the extent of + # the MALI domain, set these to None. + x_min = -630000. + x_max = 84000. + y_min = -1560000. + y_max = -860000. + + # distance from ice margin to cull (km). + # Set to a value <= 0 if you do not want + # to cull based on distance from margin. + cull_distance = 5.0 + + # mesh density parameters + # minimum cell spacing (meters) + min_spac = 1.e3 + # maximum cell spacing (meters) + max_spac = 1.e4 + # log10 of max speed (m/yr) for cell spacing + high_log_speed = 2.5 + # log10 of min speed (m/yr) for cell spacing + low_log_speed = 0.75 + # distance at which cell spacing = max_spac (meters) + high_dist = 1.e5 + # distance within which cell spacing = min_spac (meters) + low_dist = 1.e4 + + # These *_bed settings are only applied when use_bed = True. + # distance at which bed topography has no effect + high_dist_bed = 1.e5 + # distance within which bed topography has maximum effect + low_dist_bed = 5.e4 + # Bed elev beneath which cell spacing is minimized + low_bed = 50.0 + # Bed elev above which cell spacing is maximized + high_bed = 100.0 + + # mesh density functions + use_speed = True + use_dist_to_grounding_line = False + use_dist_to_edge = True + use_bed = True diff --git a/docs/developers_guide/landice/test_groups/greenland.rst b/docs/developers_guide/landice/test_groups/greenland.rst index ca464c4d86..5b27d62983 100644 --- a/docs/developers_guide/landice/test_groups/greenland.rst +++ b/docs/developers_guide/landice/test_groups/greenland.rst @@ -44,7 +44,7 @@ mesh The class :py:class:`compass.landice.tests.greenland.mesh.Mesh` defines a step for creating a variable resolution Greenland Ice Sheet mesh. -This is used by the ``high_res_mesh`` test case. +This is used by the ``mesh_gen`` test case. .. _dev_landice_greenland_smoke_test: @@ -83,9 +83,9 @@ the latter perform a 2-day restart run beginning with the end of the first. .. _dev_landice_greenland_high_res_mesh: -high_res_mesh +mesh_gen ------------- -The :py:class:`compass.landice.tests.greenland.high_res_mesh.HighResMesh` +The :py:class:`compass.landice.tests.greenland.mesh_gen.MeshGen` calls the :py:class:`compass.landice.tests.greenland.mesh.Mesh` to create the variable resolution Greenland Ice Sheet mesh. diff --git a/docs/developers_guide/landice/test_groups/thwaites.rst b/docs/developers_guide/landice/test_groups/thwaites.rst index 763e0eabc1..7217d948d4 100644 --- a/docs/developers_guide/landice/test_groups/thwaites.rst +++ b/docs/developers_guide/landice/test_groups/thwaites.rst @@ -68,11 +68,11 @@ mesh The class :py:class:`compass.landice.tests.thwaites.mesh.Mesh` defines a step for creating a variable resolution Thwaites Glacier mesh. -This is used by the ``high_res_mesh`` test case. +This is used by the ``mesh_gen`` test case. -high_res_mesh +mesh_gen ------------- -The :py:class:`compass.landice.tests.thwaites.high_res_mesh.HighResMesh` +The :py:class:`compass.landice.tests.thwaites.mesh_gen.MeshGen` calls the :py:class:`compass.landice.tests.thwaites.mesh.Mesh` to create the variable resolution Thwaites Glacier mesh. diff --git a/docs/users_guide/landice/framework.rst b/docs/users_guide/landice/framework.rst new file mode 100644 index 0000000000..89222b9a51 --- /dev/null +++ b/docs/users_guide/landice/framework.rst @@ -0,0 +1,70 @@ +.. _landice_framework: + +Land-ice Framework +================== + + +extrapolate +~~~~~~~~~~~ + +The landice framework module ``compass/landice/extrapolate.py`` provides a +function for extrapolating variables into undefined regions. It is copied +from a similar script in MPAS-Tools. + +mesh +~~~~ + +The following config options should be defined for all ``mesh_gen`` test cases (although +not necessarily with the same values shown here, which are the defaults for the 1–10km +Humboldt mesh): + +.. code-block:: cfg + + # config options for humboldt test cases + [mesh] + + # number of levels in the mesh + levels = 10 + + # Bounds of Humboldt domain. If you want the extent + # of the gridded dataset to determine the extent of + # the MALI domain, set these to None. + x_min = -630000. + x_max = 84000. + y_min = -1560000. + y_max = -860000. + + # distance from ice margin to cull (km). + # Set to a value <= 0 if you do not want + # to cull based on distance from margin. + cull_distance = 5.0 + + # mesh density parameters + # minimum cell spacing (meters) + min_spac = 1.e3 + # maximum cell spacing (meters) + max_spac = 1.e4 + # log10 of max speed (m/yr) for cell spacing + high_log_speed = 2.5 + # log10 of min speed (m/yr) for cell spacing + low_log_speed = 0.75 + # distance at which cell spacing = max_spac (meters) + high_dist = 1.e5 + # distance within which cell spacing = min_spac (meters) + low_dist = 1.e4 + + # These *_bed settings are only applied when use_bed = True. + # distance at which bed topography has no effect + high_dist_bed = 1.e5 + # distance within which bed topography has maximum effect + low_dist_bed = 5.e4 + # Bed elev beneath which cell spacing is minimized + low_bed = 50.0 + # Bed elev above which cell spacing is maximized + high_bed = 100.0 + + # mesh density functions + use_speed = True + use_dist_to_grounding_line = False + use_dist_to_edge = True + use_bed = True diff --git a/docs/users_guide/landice/framework/index.rst b/docs/users_guide/landice/framework/index.rst deleted file mode 100644 index 04e2038fae..0000000000 --- a/docs/users_guide/landice/framework/index.rst +++ /dev/null @@ -1,6 +0,0 @@ -.. _landice_framework: - -Framework -========= - -So far, the ``landice`` core has no shared framework. diff --git a/docs/users_guide/landice/index.rst b/docs/users_guide/landice/index.rst index 33efd66aec..e76ab0e4aa 100644 --- a/docs/users_guide/landice/index.rst +++ b/docs/users_guide/landice/index.rst @@ -30,4 +30,4 @@ Some helpful external links: test_groups/index suites - framework/index + framework diff --git a/docs/users_guide/landice/test_groups/greenland.rst b/docs/users_guide/landice/test_groups/greenland.rst index 7db3806fad..67188d158c 100644 --- a/docs/users_guide/landice/test_groups/greenland.rst +++ b/docs/users_guide/landice/test_groups/greenland.rst @@ -15,24 +15,24 @@ The ``landice/greenland`` test group runs tests with a coarse (20-km) The test group includes 3 test cases, each of which has one or more steps that are variants on ``run_model`` (given other names in the decomposition and restart test cases to distinguish multiple model runs), which performs time -integration of the model. There is a fourth test case, ``high_res_mesh``, that +integration of the model. There is a fourth test case, ``mesh_gen``, that creates a variable resolution Greenland Ice Sheet mesh, with the step ``mesh``. The test cases in this test group can run with either the SIA or FO velocity solvers. Running with the FO solver requires a build of MALI that includes Albany, but the SIA variant of the test can be run without Albany. The FO version uses no-slip basal boundary condition. There is no integration step -for the test case ``high_res_mesh``. +for the test case ``mesh_gen``. config options -------------- -The ``high_res_mesh`` test case uses the default config options below. +The ``mesh_gen`` test case uses the default config options below. The other test cases do not use config options. .. code-block:: cfg - [high_res_mesh] + [mesh] # number of levels in the mesh levels = 10 @@ -85,10 +85,10 @@ from a restart file saved by the first. Prognostic variables are compared between the "full" and "restart" runs at year 2 to make sure they are bit-for-bit identical. -high_res_mesh +mesh_gen ------------- -``landice/greenland/high_res_mesh`` creates a variable resolution mesh based +``landice/greenland/mesh_gen`` creates a variable resolution mesh based on the the config options listed above. This will not be the same as the pre-generated 20 km mesh used in the other three test cases because it uses a newer version of Jigsaw. Note that the basal friction optimization is diff --git a/docs/users_guide/landice/test_groups/humboldt.rst b/docs/users_guide/landice/test_groups/humboldt.rst index ad11ba8d44..40aec04a7d 100644 --- a/docs/users_guide/landice/test_groups/humboldt.rst +++ b/docs/users_guide/landice/test_groups/humboldt.rst @@ -34,7 +34,7 @@ the mesh generation options are adjusted through the config file. .. code-block:: cfg - [humboldt_mesh] + [mesh] # number of levels in the mesh levels = 10 diff --git a/docs/users_guide/landice/test_groups/kangerlussuaq.rst b/docs/users_guide/landice/test_groups/kangerlussuaq.rst index f3b31d5973..5c91bd16a6 100644 --- a/docs/users_guide/landice/test_groups/kangerlussuaq.rst +++ b/docs/users_guide/landice/test_groups/kangerlussuaq.rst @@ -24,7 +24,7 @@ the mesh generation options are adjusted through the config file. .. code-block:: cfg - [high_res_Kangerlussuaq_mesh] + [mesh] # number of levels in the mesh levels = 10 diff --git a/docs/users_guide/landice/test_groups/koge_bugt_s.rst b/docs/users_guide/landice/test_groups/koge_bugt_s.rst index cb6a6a2db6..4f65ad5009 100644 --- a/docs/users_guide/landice/test_groups/koge_bugt_s.rst +++ b/docs/users_guide/landice/test_groups/koge_bugt_s.rst @@ -24,7 +24,7 @@ the mesh generation options are adjusted through the config file. .. code-block:: cfg - [high_res_KogeBugtS_mesh] + [mesh] # number of levels in the mesh levels = 10 diff --git a/docs/users_guide/landice/test_groups/thwaites.rst b/docs/users_guide/landice/test_groups/thwaites.rst index 79b44e87e6..9a7c4bf0b1 100644 --- a/docs/users_guide/landice/test_groups/thwaites.rst +++ b/docs/users_guide/landice/test_groups/thwaites.rst @@ -27,17 +27,17 @@ integration of the model. There is a not an explicit smoke test, but the The ``decomposition_test`` and ``restart_test`` test cases in this test group can only be run with the FO velocity solvers. Running with the FO solver requires a build of MALI that includes Albany. There is no integration step for the test -case ``high_res_mesh``. +case ``mesh_gen``. config options -------------- -The ``high_res_mesh`` test case uses the default config options below. +The ``mesh_gen`` test case uses the default config options below. The other test cases do not use config options. .. code-block:: cfg - [high_res_mesh] + [mesh] # number of levels in the mesh levels = 10 @@ -84,10 +84,10 @@ second begins from a restart file saved by the first. Prognostic variables are compared between the "full" and "restart" runs to make sure they are bit-for-bit identical. -high_res_mesh +mesh_gen ------------- -``landice/thwaites/high_res_mesh`` creates a variable resolution mesh based +``landice/thwaites/mesh_gen`` creates a variable resolution mesh based on the the config options listed above. This will not be the same as the pre-generated 4-14km mesh used in ``decomposition_test`` and ``restart_test`` because it uses a newer version of Jigsaw. Note that the basal friction