.. _mesh_creation: .. |---| unicode:: U+2014 .. em dash, trimming surrounding whitespace :trim: ************* Mesh Creation ************* Uniform, Planar Meshes ====================== The most basic tool for creating an MPAS mesh is the function :py:func:`mpas_tools.planar_hex.make_planar_hex_mesh()` or the associated command-line tool ``planar_hex``. These tools create a uniform, planar mesh of hexagons that may optionally be periodic in either or both of the x and y directions. A typical call to ``make_planar_hex_mesh()`` might look like the following: .. code-block:: python from mpas_tools.planar_hex import make_planar_hex_mesh dsMesh = make_planar_hex_mesh(nx=10, ny=20, dc=1000., nonperiodic_x=False, nonperiodic_y=False) This creates a periodic mesh in both x and y that is 10 grid cells by 20 grid cells in size with 1-km resolution. The mesh is not save to a file in this example but is instead returned as an ``xarray.Dataset`` object. The command-line approach for generating the same mesh would be: .. code-block:: none $ planar_hex --nx=10 --ny=20 --dc=1000. \ --outFileName=periodic_mesh_10x20_1km.nc In this case, the resulting mesh is written to a file. The default is to make a periodic mesh (awkwardly, it's "not non-periodic"). If you want a mesh that is not periodic in one or both directions, you first need to create a mesh and then cull the cells along the periodic boundary. It doesn't hurt to run the mesh converter as well, just to be sure the mesh conforms correctly to the MPAS mesh specification: .. code-block:: python from mpas_tools.planar_hex import make_planar_hex_mesh from mpas_tools.mesh.conversion import cull, convert dsMesh = make_planar_hex_mesh(nx=10, ny=20, dc=1000., nonperiodic_x=True, nonperiodic_y=True) dsMesh = cull(dsMesh) dsMesh = convert(dsMesh) On the command line, this looks like: .. code-block:: none $ planar_hex --nx=10 --ny=20 --dc=1000. --npx --npy \ --outFileName=mesh_to_cull.nc $ MpasCellCuller.x mesh_to_cull.nc culled_mesh.nc $ MpasMeshConverter.x culled_mesh.nc nonperiodic_mesh_10x20_1km.nc Building a JIGSAW Mesh ====================== The :py:mod:`mpas_tools.mesh.creation.build_mesh` module is used create an MPAS mesh using the `JIGSAW `_ and `JIGSAW-Python (jigsawpy) `_ packages. Spherical Meshes ---------------- Spherical meshes are constructed with the function :py:func:`mpas_tools.mesh.creation.build_mesh.build_spherical_mesh()`. The user provides a 2D array ``cellWidth`` of cell sizes in kilometers along 1D arrays for the longitude and latitude (the cell widths must be on a lon/lat tensor grid) and the radius of the earth in meters. The result is an MPAS mesh file, called ``base_mesh.nc`` by default, as well as several intermediate files: ``mesh.log``, ``mesh-HFUN.msh``, ``mesh.jig``, ``mesh-MESH.msh``, ``mesh.msh``, and ``mesh_triangles.nc``. Here is a simple example script for creating a uniform MPAS mesh with 240-km resolution: .. code-block:: python #!/usr/bin/env python import numpy as np from mpas_tools.ocean import build_spherical_mesh def cellWidthVsLatLon(): """ Create cell width array for this mesh on a regular latitude-longitude grid. Returns ------- cellWidth : ndarray m x n array of cell width in km lon : ndarray longitude in degrees (length n and between -180 and 180) lat : ndarray longitude in degrees (length m and between -90 and 90) """ dlat = 10 dlon = 10 constantCellWidth = 240 nlat = int(180/dlat) + 1 nlon = int(360/dlon) + 1 lat = np.linspace(-90., 90., nlat) lon = np.linspace(-180., 180., nlon) cellWidth = constantCellWidth * np.ones((lat.size, lon.size)) return cellWidth, lon, lat def main(): cellWidth, lon, lat = cellWidthVsLatLon() build_spherical_mesh(cellWidth, lon, lat, out_filename='base_mesh.nc') if __name__ == '__main__': main() We define the resolution on a coarse (10 degree by 10 degree) grid because it is uniform. Meshes with more complex variation may require higher resolution grids to cell widths. Planar Meshes ------------- Planar meshes can be constructed with the function :py:func:`mpas_tools.mesh.creation.build_mesh.build_planar_mesh()`. Provide this function with a 2D array ``cellWidth`` of cell sizes in kilometers and 1D arrays for x and y (the cell widths must be on a 2D tensor grid). Planar meshes also require ``geom_points``, a list of point coordinates for bounding polygon for the planar mesh, and ``geom_edges``, a list of edges between points in ``geom_points`` that define the bounding polygon. As for spehrical meshes, the result is an MPAS mesh file, called ``base_mesh.nc`` by default, as well as several intermediate files: ``mesh.log``, ``mesh-HFUN.msh``, ``mesh.jig``, ``mesh-MESH.msh``, ``mesh.msh``, and ``mesh_triangles.nc``. JIGSAW Driver ------------- Underlying both spherical and planar mesh creation is the JIGSAW driver function :py:func:`mpas_tools.mesh.creation.jigsaw_driver.jigsaw_driver()`. This function is used to setup data structures and then build a JIGSAW mesh using ``jigsawpy``. Mesh Definition Tools ===================== The :py:mod:`mpas_tools.mesh.creation.mesh_definition_tools` module includes several tools for defining the ``cellWidth`` variable. Merging Cell Widths ------------------- The function :py:func:`mpas_tools.mesh.creation.mesh_definition_tools.mergeCellWidthVsLat()` is used to combine two cell-width distributions that are functions of latitude only and which asymptote to different constant values north and south of a given transition latitude with a ``tanh`` function of a given characteristic width. For example, the following code snippet will produce cell widths as a function of latitude of about 30 km south of the Arctic Circle and 10 km north of that latitude, transitioning over a characteristic "distance" of about 5 degrees. .. code-block:: python import numpy from mpas_tools.mesh.creation.mesh_definition_tools import \ mergeCellWidthVsLat lat = numpy.linspace(-90., 90., 181) cellWidthInSouth = 30. cellWidthInNorth = 10. latTransition = 66.5 latWidthTransition = 5. cellWidths = mergeCellWidthVsLat(lat, cellWidthInSouth, cellWidthInNorth, latTransition, latWidthTransition) .. _ec_mesh: Defining an Eddy-closure Mesh ----------------------------- One of the commonly used flavor of MPAS-Ocean and MPAS-Seaice meshes is designed with relatively coarse resolution in mind (requiring parameterization of ocean eddies with an "eddy closure"). This flavor of mesh has resolution that is purely a function of latitude, with 5 regions of relatively uniform resolution (north polar, northern mid-latitudes, equatorial, southern mid-latitudes and south polar) with smooth (``tanh``) transitions between these resolutions. The default EC mesh has resolutions of 35 km at the poles, 60 km at mid-latitudes and 30 km at the equator. Transitions between equatorial and mid-latitude regions are at 15 degrees N/S latitude and transitions between mid-latitude and polar regions are at 73 degrees N/S latitude. The transition near the equator is somewhat more abrupt (~6 degrees) than near the poles (~9 degrees). The switch between the transitional ``tanh`` functions is made at 40 degrees N/S latitude, where the resolution is nearly constant and no appreciable discontinuity arises. The default EC mesh can be obtained with the function :py:func:`mpas_tools.mesh.creation.mesh_definition_tools.EC_CellWidthVsLat()`: .. code-block:: python import numpy from mpas_tools.mesh.creation.mesh_definition_tools import \ EC_CellWidthVsLat lat = numpy.linspace(-90., 90., 181) cellWidths = EC_CellWidthVsLat(lat) .. _rrs_mesh: Defining a Rossby-radius Mesh ----------------------------- Another common flavor of MPAS-Ocean and MPAS-Seaice meshes is designed for higher resolutions, where the Rossby radius of deformation can be (at least partially) resolved. These meshes approximately scale their resolution in proportion to the Rossby radius. A typical Rossby Radius Scaling (RRS) mesh has a resolution at the poles that is three times finer than the resolution at the equator. For example, the RRS mesh used in E3SMv1 high resolution simulations would be defined, using the function :py:func:`mpas_tools.mesh.creation.mesh_definition_tools.RRS_CellWidthVsLat()` by: .. code-block:: python import numpy from mpas_tools.mesh.creation.mesh_definition_tools import \ RRS_CellWidthVsLat lat = numpy.linspace(-90., 90., 181) cellWidths = RRS_CellWidthVsLat(lat, cellWidthEq=18., cellWidthPole=6.) Defining an Atlantic/Pacific Mesh --------------------------------- The function :py:func:`mpas_tools.mesh.creation.mesh_definition_tools.AtlanticPacificGrid()` can be used to define a mesh that has two different, constant resolutions in the Atlantic and Pacific Oceans. Signed Distance Functions ========================= The :py:mod:`mpas_tools.mesh.creation.signed_distance` module includes several functions for creating ``cellWidth`` variables based on the signed distance from a boundary curve on the sphere. A signed distance function is positive outside the bounding shape and negative inside, with a value proportional to the distance to the nearest point on the curve (so the function is equal to zero on the curve). Signed distance functions provide a useful way ot define transitions in resolution based on complex shapes that can be defined using `geojson `_ files. These files can be created by hand, e.g. at `geojson.io `_ or in python using libraries like `shapely `_. Calls to the functions in this module require a `FeatureCollection `_ object from the `geometric_features `_ package. The ``FeatureColleciton`` must define one or more regions on the sphere from which the distance, mask, or signed distance will be computed. The ``FeatureColleciton`` could come from the predefined features included in ``geometric_features``, could be read in from a ``geojson`` file (see `Reading in Features `_), or could be created as part of a python script with ``shapely`` or other tools. In this example, we first define a base resolution using the default EC mesh (see :ref:`ec_mesh`) and then use :py:func:`mpas_tools.mesh.creation.signed_distance.signed_distance_from_geojson()` to create a signed distance function from a ``FeatureCollection`` read in from `this geojson file `_. The signed distance function is used to define a region of high resolution (12 km) around Antarctica. .. code-block:: python import numpy as np import mpas_tools.mesh.creation.mesh_definition_tools as mdt from mpas_tools.mesh.creation.signed_distance import \ signed_distance_from_geojson from geometric_features import read_feature_collection from mpas_tools.cime.constants import constants dlon = 0.1 dlat = dlon earth_radius = constants['SHR_CONST_REARTH'] nlon = int(360./dlon) + 1 nlat = int(180./dlat) + 1 lon = np.linspace(-180., 180., nlon) lat = np.linspace(-90., 90., nlat) cellWidth = mdt.EC_CellWidthVsLat(lat) # broadcast cellWidth to 2D _, cellWidth = np.meshgrid(lon, cellWidthVsLat) # now, add the high-res region fc = read_feature_collection('high_res_region.geojson') so_signed_distance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) # Equivalent to 20 degrees latitude trans_width = 1600e3 trans_start = -500e3 dx_min = 12. weights = 0.5 * (1 + np.tanh((so_signed_distance - trans_start) / trans_width)) cellWidth = dx_min * (1 - weights) + cellWidth * weights Sometimes it can be useful to extract just the mask of the region of interest (defined as ``0`` outside the the region and ``1`` inside it) or the unsigned distance. For these purposes, use the functions :py:func:`mpas_tools.mesh.creation.signed_distance.mask_from_geojson()` and :py:func:`mpas_tools.mesh.creation.signed_distance.distance_from_geojson()`, respectively. The grid does not necessarily have to be uniform in resolution. Here is a relatively complicated example where the highest resolution (1/100 degree) of the lon/lat grid is concentrated within +/- 10 degrees of 0 degrees longitude and 0 degrees latitude and the resolution is ~2 degrees elsewhere. The shape ``test.geojson`` is given below, and is in the high-res region: .. code-block:: python import numpy import matplotlib.pyplot as plt from geometric_features import read_feature_collection from mpas_tools.cime.constants import constants from mpas_tools.mesh.creation.signed_distance import \ distance_from_geojson, mask_from_geojson, signed_distance_from_geojson def gaussian(dcoarse, dfine, cmin, cmax, center, width, iter_count=10): def get_res(xval): res = dcoarse + \ (dfine - dcoarse) * numpy.exp(-(xval - center) ** 2 / (2. * width**2)) return res x = [cmin] while x[-1] < cmax: d = get_res(x[-1]) for index in range(iter_count): x_mid = x[-1] + 0.5*d d = get_res(x_mid) x.append(x[-1] + d) x = numpy.array(x) factor = (cmax - cmin)/(x[-1] - x[0]) x = cmin + factor*(x - x[0]) return x def main(): dcoarse = 2. dfine = 0.01 width = 10. lon = gaussian(dcoarse=dcoarse, dfine=dfine, cmin=-180., cmax=180., center=0., width=width) lat = gaussian(dcoarse=dcoarse, dfine=dfine, cmin=-90., cmax=90., center=0., width=width) earth_radius = constants['SHR_CONST_REARTH'] fc = read_feature_collection('test.geojson') distance = distance_from_geojson(fc, lon, lat, earth_radius, max_length=dfine) mask = mask_from_geojson(fc, lon, lat, max_length=dfine) signed_distance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=dfine) Here is ``test.geojson``: .. code-block:: python { "type": "FeatureCollection", "features": [ { "type": "Feature", "properties": { "name": "weird shape", "component": "ocean", "object": "region", "author": "Xylar Asay-Davis" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 3.8232421875000173, 3.283113991917228 ], [ 1.5600585937500142, 3.8752158957336085 ], [ 0.9448242187500171, 1.8453839885731744 ], [ -0.3186035156249841, 0.8239462091017558 ], [ 0.6372070312500188, -0.37353251022881745 ], [ 2.3181152343750147, 0.03295898255728466 ], [ 3.636474609375017, -0.3405741662837591 ], [ 4.163818359375015, 1.5159363834516735 ], [ 2.9443359375000164, 1.9771465537125645 ], [ 3.8232421875000173, 3.283113991917228 ] ] ] } } ] }