Mesh Creation

Uniform, Planar Meshes

The most basic tool for creating an MPAS mesh is the function 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:

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:

$ 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:

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:

$ 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 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 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:

#!/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 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 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 mpas_tools.mesh.creation.mesh_definition_tools module includes several tools for defining the cellWidth variable.

Merging Cell Widths

The function 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.

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)

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 mpas_tools.mesh.creation.mesh_definition_tools.EC_CellWidthVsLat():

import numpy
from mpas_tools.mesh.creation.mesh_definition_tools import \
    EC_CellWidthVsLat

lat = numpy.linspace(-90., 90., 181)
cellWidths = EC_CellWidthVsLat(lat)

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 mpas_tools.mesh.creation.mesh_definition_tools.RRS_CellWidthVsLat() by:

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 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 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 Defining an Eddy-closure Mesh) and then use 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.

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 mpas_tools.mesh.creation.signed_distance.mask_from_geojson() and 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:

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:

{
  "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
            ]
          ]
        ]
      }
    }
  ]
}