Source code for compass.landice.tests.humboldt.mesh

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 import write_netcdf
from mpas_tools.logging import check_call

from compass.step import Step
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

[docs]class Mesh(Step): """ A step for creating a mesh and initial condition for humboldt test cases Attributes ---------- mesh_type : str The resolution or mesh type of the test case """
[docs] def __init__(self, test_case): """ Create the step Parameters ---------- test_case : compass.TestCase The test case this step belongs to mesh_type : str The resolution or mesh type of the test case """ super().__init__(test_case=test_case, name='mesh') self.add_output_file(filename='') self.add_output_file(filename='') self.add_input_file( filename='', target='', database='') self.add_input_file(filename='Humboldt.geojson', package='compass.landice.tests.humboldt', target='Humboldt.geojson', database=None) self.add_input_file(filename='', target='', database='')
# no setup() method is needed
[docs] def run(self): """ Run this step of the test case """ logger = self.logger config = self.config section = config['humboldt']'calling build_cell_wdith') cell_width, x1, y1, geom_points, geom_edges = self.build_cell_width()'calling build_planar_mesh') build_planar_mesh(cell_width, x1, y1, geom_points, geom_edges, logger=logger) dsMesh = xarray.open_dataset('')'culling mesh') dsMesh = cull(dsMesh, logger=logger)'converting to MPAS mesh') dsMesh = convert(dsMesh, logger=logger)'writing') write_netcdf(dsMesh, '') # If no number of levels specified in config file, use 10 levels = section.get('levels')'calling') args = ['', '-i', '', '-o', '', '-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.'calling') args = ['', '-s', '', '-d', '', '-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.:'calling') args = ['', '-f', '', '-m' 'distance', '-d', cullDistance] check_call(args, logger=logger) else:'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.'calling') args = ['', '-f', '', '-p', 'gis-gimp'] check_call(args, logger=logger)'calling MpasMaskCreator.x') args = ['MpasMaskCreator.x', '', '', '-f', 'Humboldt.geojson'] check_call(args, logger=logger)'culling to geojson file') dsMesh = xarray.open_dataset('') humboldtMask = xarray.open_dataset('') dsMesh = cull(dsMesh, dsInverse=humboldtMask, logger=logger) write_netcdf(dsMesh, '')'Marking horns for culling') args = ['', '-f', ''] check_call(args, logger=logger)'culling and converting') dsMesh = xarray.open_dataset('') dsMesh = cull(dsMesh, logger=logger) dsMesh = convert(dsMesh, logger=logger) write_netcdf(dsMesh, '')'calling') args = ['', '-i', '', '-o', '', '-l', levels, '-v', 'glimmer', '--beta', '--thermal', '--obs', '--diri'] check_call(args, logger=logger)'calling') args = ['', '-s', '', '-d', '', '-m', 'b', '-t'] check_call(args, logger=logger)'Marking domain boundaries dirichlet') args = ['', '-f', ''] check_call(args, logger=logger)'calling') args = ['', '-f', '', '-p', 'gis-gimp'] check_call(args, logger=logger)'creating') make_graph_file(mesh_filename='', graph_filename='')
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('', '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, window_size=1.e5) # optional - plot distance calculation # plt.pcolor(distToEdge/1000.0); plt.colorbar(); # Set cell widths based on mesh parameters set in config file cell_width = set_cell_width(self, section='humboldt', thk=thk, vx=vx, vy=vy, dist_to_edge=distToEdge, dist_to_grounding_line=None) # plt.pcolor(cell_width); plt.colorbar(); return (cell_width.astype('float64'), x1.astype('float64'), y1.astype('float64'), geom_points, geom_edges)