Skip to content

Conversation

@alexolinhager
Copy link

This PR creates the option (max_res_in_ocn=True) to maintain the highest resolution used on land ice cells in all ocean cells when generating the GIS MALI mesh. This is done to resolve fjord geometry and coastal bathymetry necessary to translate ocean properties to glacier termini when coupling MALI to MPAS-Ocean (or prescribing ocean forcing to MALI in standalone mode). As this has the potential to drastically increase mesh generation time, so there is an additional option define_bnds_by_geojson=True that performs an initial cull of the BedMachine dataset according to the bounds set by compass/landice/tests/greenland/gis_contShelfExtent_PS.geojson before generating the mesh in Jigsaw.

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

high_bed = section.getfloat('high_bed')

# convert km to m
cull_distance = section.getfloat('cull_distance') * 1.e3
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We definitely don't want to remove culling by distance entirely.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This call of cull_distance is only used to assign large widths for cells that are to be culled later on (to speed up efficiency). I think Matt and I decided this created more problems than it was actually worth. The actual culling still happens here:
https://github.com/MPAS-Dev/compass/blob/b3d1faba9502e38209b3376b2241c35b9dbffa54/compass/landice/mesh.py#L713C5-L722C63

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, thanks for pointing that out. I don't know that I necessarily agree with this choice as a general decision for all meshes going forward. I can imagine a ton of wasted compute time trying to make a high-resolution AIS or GrIS mesh that creates high-resolution cells hundreds of km away from the ice sheet, only to cull them away afterwards. I think we need some demonstration that this doesn't increase cost in general too much.

@alexolinhager alexolinhager requested a review from trhille October 27, 2025 20:26
@trhille
Copy link
Collaborator

trhille commented Oct 27, 2025

@alexolinhager, this is missing an updated mesh_gen.cfg file.

Copy link
Collaborator

@trhille trhille left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@alexolinhager, thanks for this really useful PR. I have a number of questions and thoughts, mostly regarding generalization. I've also pointed out a few logical inconsistencies that should be easy to fix.

# thermal forcing parameterizations. Otherwize, telescope out
# to coarse resolution with distance from the ice front.
if section.get('max_res_in_ocn') == 'True':
spacing_bed[np.logical_and(bed < 0, thk == 0)] = min_spac
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems restrictive to impose that the ocean resolution is the same as the minimum ice-sheet resolution. If the ocean model resolution is >10km and we have ≤1 km minimum ice-sheet resolution, then that seems like a lot of wasted cells. What about explicitly setting the desired resolution in the open ocean and then either defining a length-scale over which the resolution varies from the ice edge outwards, or defining an open-ocean mask with a geojson file?

spacing_bed[np.logical_and(bed < 0, thk == 0)] = min_spac
print("Maximizing resolution in ocean. {} ocean cells \
".format(np.sum(np.logical_and(bed < 0, thk == 0))))
else:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should not be within an else-block, or we will not be able to set maximum resolution in the ocean and have bed topography influence resolution for the ice sheet.

high_bed = section.getfloat('high_bed')

# convert km to m
cull_distance = section.getfloat('cull_distance') * 1.e3
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, thanks for pointing that out. I don't know that I necessarily agree with this choice as a general decision for all meshes going forward. I can imagine a ton of wasted compute time trying to make a high-resolution AIS or GrIS mesh that creates high-resolution cells hundreds of km away from the ice sheet, only to cull them away afterwards. I think we need some demonstration that this doesn't increase cost in general too much.

bnd = section.get(option)
if bnd != 'None':
bnds[index] = float(bnd)
# If necessary, get bounds defined by user or use bound of gridded dataset
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm confused about why this needs to be changed.


return (cell_width.astype('float64'), x1.astype('float64'),
y1.astype('float64'), geom_points, geom_edges, flood_mask)
if not calc_geom_bnds:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's a bad idea to return a different size tuple based on the input arguments. If it's really necessary to put the if-statement above, it would be better to have an else-statement that just sets geom_points, geom_edges = None, None and keep the return statement simple.

self, section_name=section_name,
gridded_dataset=source_gridded_dataset_2km,
flood_fill_start=[100, 700])
if section_gis.get("define_bnds_by_geojson") == 'True':
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As noted above, I'd rather that build_cell_width return None, None for geom_points, geom_edges here, and then move the call set_geojson_geom_points_and_edges() after the call to build_cell_width to get those arrays.

ds.to_netcdf(self.mesh_filename, 'a')


def set_geojson_geom_points_and_edges(self):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be generalized and moved to the landice.mesh module.

andrewdnolan and others added 7 commits October 28, 2025 14:46
Creates the option to maintain resolution in the ocean. Used to resolve
fjords and bathymetry when coupling to MPAS-Ocean.
This commit creates the option to maintain high resolution in all ocean
cells (same resolution as what is prescribed at outlet glacier termini).
This is necessary to resolve the fjord geometry that thermal forcing
parameterizations rely on. To reduce time when generating these higher
resolution meshes, this commit also introduces an option to do an
initial cull of the bedMachine dataset using the extent of the
contintental shelf before generating the GIS mesh.
This commit moves the calculation of geom_points and geom_edges
using a geojson of the gis contintental shelf extent
to compass/landice/tests/greenland/mesh.py, instead of
compass/landice/mesh.py. If this options is not being used, thn
geom_points and geom_edges are calculated within build_cell_width, as
per usual
Co-authored-by: Trevor Hillebrand <trhille@lanl.gov>
Re-introduces code that sets a large cell width for cells that will be
culled by distance from ice edge. Only applied when not using a geojson
to pre-cull domain
Ensures the bedTopography constraint on mesh resolution is enforced
regardless 'max_res_in_ocn'. Instead, grid resolution in the ocean
is maximized following the bedTopgraphy constraint, overwritting it
where necessary.
@alexolinhager alexolinhager force-pushed the alexolinhager/hiRes_ocean_GIS_mesh branch from 5f0591d to 688b86e Compare October 28, 2025 21:47
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants