Create test for extracting a MALI subdomain from a larger domain

[matthewhoffman]

This merge adds a test case for a subdomain extractor. It allows a regional domain to be quickly created from an existing whole-ice-sheet domain using a mask. The regional mask can come from a region_mask MPAS netCDF file or a geojson file. Optionally, ancillary files (e.g. forcing or parameter files) can be remapped at the same time.

Checklist

• User's Guide has been updated
• Developer's Guide has been updated
• API documentation in the Developer's Guide (api.rst) has any new or modified class, method and/or functions listed
• Documentation has been built locally and changes look as expected
• Document (in a comment titled Testing in this PR) any testing that was used to verify the changes

[matthewhoffman]

Testing

I created an Amery subdomain with the following cfg file on perlmutter:

[subdomain]

# path to file from which to extract subdomain
source_file = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/AIS_4to20km_r01_20220907_m5_drop_bed_20m_bulldoze_troughs_75_to_400m_Enderby_maxstiffness_0.8_TG_pinning_40maf_bedmap2_surface_ASE_05perc_seafloor_mu.nc

# path to region mask file for source_file
region_mask_file = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/AIS_4to20km_r01_20220907.regionMask_ismip6.nc

# region number to extract
region_number = 3

# filename for the subdomain to be generated
dest_file_name = Amery.nc

# mesh projection to be used for setting lat/lon values
# Should match the projection used in the source_file
# Likely one of 'ais-bedmap2' or 'gis-gimp'
mesh_projection = ais-bedmap2

# whether to extend mesh into the open ocean along the ice-shelf margin
# This is necessary if the region mask ends right along the ice-shelf edge,
# or if the ice-shelf is covered by two regions.
# It is recommended to try extracting a subdomain with this False, and if the
# ocean buffer is inadequate, trying again with this True.
extend_ocean_buffer = True

# number of iterations to grow open ocean buffer along ice-shelf margin
# Only used if extend_ocean_buffer is True
# Should be equal to approximately the number of ocean buffer cells in the
# source_file
grow_iters = 15

Result:

Result setting extend_ocean_buffer = False:

[matthewhoffman]

Temporary docs are here: https://portal.nersc.gov/cfs/fanssie/hoffman2/compass-docs-subdomain-extractor/html/index.html
Specifically:

• https://portal.nersc.gov/cfs/fanssie/hoffman2/compass-docs-subdomain-extractor/html/users_guide/landice/test_groups/mesh_modifications.html
• https://portal.nersc.gov/cfs/fanssie/hoffman2/compass-docs-subdomain-extractor/html/developers_guide/landice/test_groups/mesh_modifications.html
• https://portal.nersc.gov/cfs/fanssie/hoffman2/compass-docs-subdomain-extractor/html/developers_guide/landice/api.html#mesh-modifications

[trhille]

@matthewhoffman, is this expected behavior for extend_mesh = False (left) vs True (right)? i.e., is it always desirable to extend only into floating ice and exclude ice rises, etc?

[trhille]

@matthewhoffman, thanks for creating this great tool! See my comments regarding a small amount of repeated code for the flood-fill, and also my question about the desired behavior of extend_mesh. I'll wait to approve until you respond to those.

[matthewhoffman]

@trhille , I addressed the comments from your review and added a few significant new features:

• ability to use geojson file instead of existing region_mask file
• ability to also remap ancillary files on the same mesh (e.g. forcing, parameter files)
• the second item required adding a new remap method using ncremap instead of our standard MALI interpolation tools
• I also moved the flood fill to landice framework

On Perlmutter this cfg can be used to test this PR out. As set up, setting region_definition = region_mask_file will generate a subdomain mesh and 3 forcing files for Amery. Changing that to region_definition = geojson will instead create a subdomain for our Thwaites Glacier 'minimal' region definition.

[subdomain]

# path to file from which to extract subdomain
source_file = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/AIS_4to20km_r01_20220907_m5_drop_bed_20m_bulldoze_troughs_75_to_400m_Enderby_maxstiffness_0.8_TG_pinning_40maf_bedmap2_surface_ASE_05perc_seafloor_mu.nc

# method for defining region
# one of 'region_mask_file' or 'geojson'
region_definition = geojson
#region_definition = region_mask_file

# path to geojson file to be used if region_definition='geojson'
geojson_file = /global/cfs/cdirs/fanssie/users/hoffman2/compass/subdomain_extractor_cleanup/compass/landice/tests/thwaites/thwaites_minimal.geojson


# path to region mask file for source_file
region_mask_file = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/AIS_4to20km_r01_20220907.regionMask_ismip6.nc

# region number to extract
region_number = 3

# filename for the subdomain to be generated
dest_file_name = Subdomain.nc

# mesh projection to be used for setting lat/lon values
# Should match the projection used in the source_file
# Likely one of 'ais-bedmap2' or 'gis-gimp'
mesh_projection = ais-bedmap2

# whether to extend mesh into the open ocean along the ice-shelf margin
# This is necessary if the region mask ends right along the ice-shelf edge,
# or if the ice-shelf is covered by two regions.
# It is recommended to try extracting a subdomain with this False, and if the
# ocean buffer is inadequate, trying again with this True.
extend_ocean_buffer = True

# number of iterations to grow open ocean buffer along ice-shelf margin
# Only used if extend_ocean_buffer is True
# Should be equal to approximately the number of ocean buffer cells in the
# source_file
grow_iters = 15

# method for performing interpolation
# 'ncremap' uses pyremap to call ESMF_RegridWeightGen to generate a
# nstd weight file and then uses ncremap to perform remapping.
# This method supports interpolating ancillary files (below)
# but likely needs to be run on a compute node and is more fragile.
# 'mali_interp' uses the MALI interpolation script interpolate_to_mpasli_grid.py
# This method does not support ancillary files but may be more robust
# and can likely be run quickly on a login node.
interp_method = ncremap

# optional forcing files that could also be interpolated
# interpolating these files requires using the 'ncremap' interp_method
extra_file1 = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/forcing/ais_mesh_4to20km_res/tier1_experiments/expAE01/AIS_4to20km_r01_20220907_smb_NorESM1-M_RCP26-repeat_2300_minus1_bare_land.nc
extra_file2 = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/forcing/ais_mesh_4to20km_res/tier1_experiments/ctrlAE/AIS_4to20km_r01_20220907_obs_TF_1995-2017_8km_x_60m.nc
extra_file3 = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/basin_and_coeff_gamma0_DeltaT_quadratic_non_local_after_relaxation.nc
extra_file4 = None
extra_file5 = None

[trhille]

Nice work, @matthewhoffman! I'm approving based on inspection and your testing. Someone else should test independently, as we discussed.

[alexolinhager]

Looked through the code and ran the compass test case on Perlmutter. Everything seems to be working as intended

[xylar]

@matthewhoffman, this should have been named dev_mesh_modifications. It has the same name as the section in the user's guide, which is producing warnings when building the documentation. Could you open a small PR to fix this when you have time?

(I accidentally commented previously on the commit, not the PR.)

[matthewhoffman]

Addressed in #808