Add functions for culling an MPAS dataset

[xylar]

This merge adds tools for culling an MPAS dataset, given a dataset on the unculled mesh as well as an unculled and culled horizontal MPAS mesh.

The merge also provides functions for making (and storing) the index mapping between the culled and unculled meshes so it can be reused.

[xylar]

Testing

I used these functions to crop an initial condition from the Compass ISOMIP+ Ocean0 test case (a problem small enough for my laptop) to only the portion of the domain without ice-shelf cavities. The code I used has been added to the documentation here as an example.

I have not tried to perform a forward run with the resulting initial condition. I will do so on Chrysalis using a global ocean configuration and will report back here soon.

[cbegeman]

@xylar This looks great! Do you have a compass branch handy that I could use for testing?

[xylar]

Do you have a compass branch handy that I could use for testing?

I am not testing this in compass, I'm just making my own local python scripts and then wiring up the symlinks manually. Here is what I did in the Icos240 case within the initial_state step:

#!/usr/bin/env python
import os

import xarray as xr

from mpas_tools.io import write_netcdf
from mpas_tools.mesh.conversion import cull
from mpas_tools.mesh.cull import write_map_culled_to_base, write_culled_dataset
from mpas_tools.logging import LoggingContext

base_mesh_filename = 'mesh.nc'
culled_mesh_filename = 'mesh_no_isc.nc'
map_filename = 'no_isc_to_culled_map.nc'

if not os.path.exists(culled_mesh_filename):
    ds_culled_mesh = xr.open_dataset(base_mesh_filename)
    ds_init = xr.open_dataset('initial_state.nc')
    ds_culled_mesh['cullCell'] = ds_init.landIceMask
    ds_culled_mesh_no_isc = cull(ds_culled_mesh,
                                 graphInfoFileName='mesh_no_isc_graph.info')
    write_netcdf(ds_culled_mesh_no_isc, culled_mesh_filename)

if not os.path.exists(map_filename):
    write_map_culled_to_base(base_mesh_filename=base_mesh_filename,
                             culled_mesh_filename=culled_mesh_filename,
                             out_filename=map_filename)

with LoggingContext('test') as logger:
    for prefix in ['initial_state', 'init_mode_forcing_data']:

        in_filename = f'{prefix}.nc'
        out_filename = f'{prefix}_no_isc.nc'

        write_culled_dataset(in_filename=in_filename, out_filename=out_filename,
                             base_mesh_filename=base_mesh_filename,
                             culled_mesh_filename=culled_mesh_filename,
                             map_culled_to_base_filename=map_filename,
                             logger=logger)

Then, I manually modified the symlinks in the init/ssh_adjustment and performance_test/prognostic_ice_shelf_melt steps to point to the culled files.

The process for the Icos30 restart run is even more manual, since I am continuing that run from the dynamic-adjustment simulation I used to make the IcoswISC30E3r5 mesh in the first place (to ensure that the mesh is exactly the right one). I can describe the process in more detail but I envisioned it just as a clunky one-off proof of concept so I'm not sure how worthwhile that will be.

[xylar]

Here is the script I used to cull an IcoswISC30E3r5 restart file:

#!/usr/bin/env python
import os

import xarray as xr

from mpas_tools.io import write_netcdf
from mpas_tools.mesh.conversion import cull
from mpas_tools.mesh.cull import write_map_culled_to_base, write_culled_dataset
from mpas_tools.logging import LoggingContext, check_call

with LoggingContext('test') as logger:

    restart_filename = '../unculled/rst.0001-01-21_00.00.00.nc'
    base_mesh_filename = restart_filename
    culled_mesh_filename = 'mesh_no_isc.nc'
    map_filename = 'no_isc_to_culled_map.nc'

    if not os.path.exists(culled_mesh_filename):
        mesh_vars = [
            'areaCell', 'cellsOnCell', 'edgesOnCell', 'indexToCellID',
            'latCell', 'lonCell', 'meshDensity', 'nEdgesOnCell',
            'verticesOnCell', 'xCell', 'yCell', 'zCell', 'angleEdge',
            'cellsOnEdge', 'dcEdge', 'dvEdge', 'edgesOnEdge',
            'indexToEdgeID', 'latEdge', 'lonEdge', 'nEdgesOnCell',
            'nEdgesOnEdge', 'verticesOnEdge', 'weightsOnEdge', 'xEdge',
            'yEdge', 'zEdge', 'areaTriangle', 'cellsOnVertex', 'edgesOnVertex',
            'indexToVertexID', 'kiteAreasOnVertex', 'latVertex',
            'lonVertex', 'xVertex', 'yVertex', 'zVertex']
        ds_unculled_mesh = xr.open_dataset(base_mesh_filename)
        ds_unculled_mesh = ds_unculled_mesh[mesh_vars]
        ds_init = xr.open_dataset(restart_filename)
        ds_unculled_mesh['cullCell'] = ds_init.landIceMask

        write_netcdf(ds_unculled_mesh, 'unculled_mesh.nc')

        args = ['MpasCellCuller.x', 'unculled_mesh.nc', culled_mesh_filename]
        check_call(args=args, logger=logger)

    if not os.path.exists(map_filename):
        write_map_culled_to_base(base_mesh_filename=base_mesh_filename,
                                 culled_mesh_filename=culled_mesh_filename,
                                 out_filename=map_filename)

    in_filename = restart_filename
    out_filename = 'rst.0001-01-21_00.00.00.nc'

    write_culled_dataset(in_filename=in_filename, out_filename=out_filename,
                         base_mesh_filename=base_mesh_filename,
                         culled_mesh_filename=culled_mesh_filename,
                         map_culled_to_base_filename=map_filename,
                         logger=logger)

    in_filename = '../unculled/init_mode_forcing_data.nc'
    out_filename = 'init_mode_forcing_data_no_isc.nc'

    write_culled_dataset(in_filename=in_filename, out_filename=out_filename,
                         base_mesh_filename=base_mesh_filename,
                         culled_mesh_filename=culled_mesh_filename,
                         map_culled_to_base_filename=map_filename,
                         logger=logger)

I made the following symlinks in an unculled directory:

init_mode_forcing_data.nc -> /lcrc/group/e3sm/ac.xylar/compass_1.2/chrysalis/e3smv3-meshes/icoswisc30e3r5/ocean/global_ocean/IcoswISC/WOA23/init/initial_state/init_mode_forcing_data.nc
rst.0001-01-21_00.00.00.nc -> /lcrc/group/e3sm/ac.xylar/compass_1.2/chrysalis/e3smv3-meshes/icoswisc30e3r5/ocean/global_ocean/IcoswISC/WOA23/dynamic_adjustment/restarts/rst.0001-01-21_00.00.00.nc

I didn't have the original culled mesh (the one that still has ice-shelf cavities but not land) so I had to extract it from the restart file.

[cbegeman]

@xylar Great. Thanks for sharing the scripts as well. By visual inspection and looking at the netcdf files in your directories, this looks good and makes sense to me. Would you like me to do any additional testing?

[xylar]

@cbegeman, I think my testing is sufficient. Thanks for having a look. You could also take a look at the resulting in:

/lcrc/group/e3sm/ac.xylar/compass_1.2/chrysalis/e3smv3-meshes/icosxisc30e3r7/ocean/utility/cull_restarts/cull

and do any sanity checks that come to mind. I haven't yet looked at them in ParaView, for example and that would very likely be a good idea.

[xylar]

Thanks @cbegeman!