Add python code for computing planar hex meshes

[xylar]

Currently produces bit-for-bit identical results to periodic_hex/periodic_grid but could be extended in the future to produce non-periodic meshes.

The reason for this update is that the python code will likely be much easier to incorporate into workflows like COMPASS because it doesn't require compilation or the Fortran NetCDF library. It could also be incorporated quite easily into a conda package for mesh tools (see #243 and #245).

partially addresses #245

[xylar]

Testing

I tested both a 10x20 cell, 1km mesh and a 500x500 cell, 10km mesh. Both were bit-for-bit identical to results from periodic_hex, tested using the function here:

MPAS-Tools/mesh_tools/planar_hex/planar_hex.py

Lines 349 to 372 in f1ec117

	def make_diff(mesh, refMeshFileName, diffFileName):
	refMesh = xarray.open_dataset(refMeshFileName)
	diff = xarray.Dataset()
	for variable in mesh.data_vars:
	if variable in refMesh:
	diff[variable] = mesh[variable] - refMesh[variable]
	print(diff[variable].name, float(numpy.abs(diff[variable]).max()))
	else:
	print('mesh has extra variable {}'.format(mesh[variable].name))
	for variable in refMesh.data_vars:
	if variable not in mesh:
	print('mesh mising variable {}'.format(refMesh[variable].name))
	for attr in refMesh.attrs:
	if attr not in mesh.attrs:
	print('mesh mising attribute {}'.format(attr))
	for attr in mesh.attrs:
	if attr not in refMesh.attrs:
	print('mesh has extra attribute {}'.format(attr))
	write_netcdf(diff, diffFileName)

Surprisingly, the python code is 8x faster for the larger mesh than the Fortran code.

[xylar]

I chose not to have the python code add the Time dimension or any tracers or vertical levels to the mesh file. It looks like periodic_hex also stopped adding tracers and vertical levels but namelist.input still has options for them. I don't believe the empty Time dimension is useful since no variables have that dimension.

I also chose not to produce the graph.info.part.* files. My guess is that any workflow we have typically culls cells and then produces its own partition files. But I would like to know if anyone uses that functionality and if I should add it back in.

[matthewhoffman]

@xylar , I confirmed this works as advertised, and I agree with the decision to leave out Time and any non-mesh variables. I also agree that there is no real value in the behavior of periodic_hex to output graph.info.part.* files, and the new script does not need to do that. (The only really useful file is graph.info, but periodic_hex currently does not output that, and there is an existing script in the MPAS-Tools repo that can be used to create it.) I also confirmed that this is substantially faster than the current code!

A couple of things I noticed:

1. Is requiring the xarray module necessary? It seems like all functionality could be produced using just the netCDF4 module (which you already have as a requirement).

2. In periodic_hex, we added the requirement that ny be an even number, otherwise you get a mesh where the y-periodicity is messed up. That should be added here too.

[xylar]

@matthewhoffman, thanks for the review!

1. Is requiring the xarray module necessary? It seems like all functionality could be produced using just the netCDF4 module

Yes, I rely extensively on xarray to handle the data structures and math between variables here. So removing it would require a substantial rewrite, which I don't think is justified. xarray is very easy to install with conda. A lot of the speed-up likely comes because of xarray.

2. In periodic_hex, we added the requirement that ny be an even number, otherwise you get a mesh where the y-periodicity is messed up. That should be added here too.

I must have missed that check. I agree that it's necessary.

[xylar]

I had a discussion with @matthewhoffman. Both of us would like to know how others feel about the dependency on xarray.

[milenaveneziani]

For me, it's a no brainer: if it works, and it's fast, and removing xarray means a substantial rewriting of the code, then why remove it??

[xylar]

Thanks or weighing in @milenaveneziani. @mgduda, do you have an opinion about this? Keep in mind that my primary reason for developing it is as a tool that can more easily be incorporated into a conda package (which can include xarray as a dependency without any extra pain for the user). This need not be a replacement for periodic_hex, just an equivalent tool.

[xylar]

@maltrud had raised the concern that xarray output wasn't working well in MPAS for him and @bradyrx. I can confirm that I can run the results of this tool through MpasMeshConverter.x so it should at least be possible to start with this tool, use the mesh converter to "fix" it for MPAS if necessary, and proceed from there. I haven't verified either way if I can run MPAS with the meshes this tool produces. That'd definitely be a prerequisite of merging this PR.

[bradyrx]

@xylar, to clarify, the issue I had was with using xarray to modify particle initialization files in LIGHT. For instance, using xarray to quickly filter out particles north of some latitude and then saving with .to_netcdf() would cause netCDF (I/O) errors when running MPAS.

Even taking care to set the output format to NetCDF3_64BIT wouldn't help. I.e., metadata of the modified netCDF looked identical to that of files that worked with MPAS, but something under the hood from xarray caused bugs. Note that I did not try MpasMeshConverter.x to fix the bug. We went with the netCDF4 package in python to fix the issues.

[matthewhoffman]

@xylar , there was no assumption before that output from perdiodic_hex should work as a valid MPAS mesh without running through the mesh converter tool first. So I think it's perfectly fine (and preferable) to still require that.

[xylar]

@bradyrx, your case wouldn't have worked with the mesh converter, I suspect. It's only for MPAS meshes, not any kind of MPAS data (e.g. particles). I have had (minor) issues with fill values from xarray in my other work, which is why I have a helper function to use standard fill values in the output. It just occurs to me that I have a run that's mysteriously crashing, and xarray could be the cause! I wonder what it does that MPAS doesn't like. Definitely worth looking into but it's outside the scope of this PR, for sure.

[xylar]

there was no assumption before that output from perdiodic_hex should work as a valid MPAS mesh without running through the mesh converter tool first. So I think it's perfectly fine (and preferable) to still require that.

@matthewhoffman, yep, that sound reasonable. If results from periodic_hex didn't work without running mesh converter, there's no reason this tool would (or should) be any different.

[matthewhoffman]

To clarify, the thinking Doug always had was to make only one tool (mesh converter) responsible for precise implementation of the MPAS mesh spec. In the early days we kept running into issues where such and such tool or MPAS framework routine would appear broken only to discover some obscure detail of the mesh spec was not being obeyed in a mesh. The mesh spec is complicated enough that requiring every tool to follow it to the letter isn't worth the effort.

[xylar]

That's a good guiding principle, for sure. It's obviously a little problematic that the tool for enforcing the mesh spec was written using c++, a language used by very few of the existing developers, and using a version of the NetCDF c++ library that was short-lived. But hindsight is always 20/20 and we have ways of building that tool wherever needed now (#243). For the record, I tried rewriting the mesh converter in python and it was a nightmare (and very slow) so I have a lot of respect of Doug's work on that tool.

[mark-petersen]

I tested with this COMPASS test case:

ocean/baroclinic_channel/10km/restart_test

Except I substituted this for init_step1:

./planar_hex.py --nx 16 --ny 50 --dc 10000 -o base_mesh.nc

It worked great! The other init and forward steps went all the way through based on this mesh, with no problem.

[mark-petersen]

@siddharthabishnu @amrapallig and @ipdemes This is a very useful tool that Xylar just added to create a 3D regular hex domain with your specifications. Please give it a try. It is in this branch of MPAS-Tools right now, but will be merged soon. With this tool, you can do a spatial convergence test.

Here is a quick way to access it while it is on this branch:

git clone git@github.com:xylar/MPAS-Tools.git xylar-MPAS-Tools
cd xylar-MPAS-Tools
git checkout -b python_planar_hex origin/python_planar_hex
cd mesh_tools/planar_hex/
ls -l
-rwxr-xr-x 1 mpeterse mpeterse 16087 Feb 28 08:52 planar_hex.py

then use it with:

./planar_hex.py --nx 16 --ny 50 --dc 10000 -o base_mesh.nc

You will need the right python environment. On LANL IC:

module unload python; source /usr/projects/climate/SHARED_CLIMATE/anaconda_envs/load_latest_e3sm_unified.sh

If you have comments, or problems, you can post on this PR.

[mark-petersen]

@siddharthabishnu @amrapallig and @ipdemes I should also say, after this command:

./planar_hex.py --nx 16 --ny 50 --dc 10000 -o base_mesh.nc

you still need to run the c-code mpas mesh converter:

./MpasMeshConverter.x base_mesh.nc mesh.nc

which produces the graph.info file. That code is housed here: https://github.com/MPAS-Dev/MPAS-Tools/tree/master/mesh_tools/mesh_conversion_tools.

@xylar Does planar_hex.py produce a valid MPAS mesh, or does MpasMeshConverter always need to be run as a second step?

[xylar]

@mark-petersen, see discussion above. @matthewhoffman thinks that no tool other than MpasMeshConverter.x should be trusted to produce a valid MPAS mesh, so we should always run any mesh through that tool. I can't speak for whether things work without that for planar_hex.py but I certainly wouldn't guarantee it.

[mgduda]

@matthewhoffman Actually, we do assume that the output of periodic_hex can be directly used by MPAS-Atmosphere for idealized simulations, with no need to run through a mesh converter tool. I can definitely help in testing that the output of the proposed Python tool works in this regard.

[xylar]

@mgduda

I can definitely help in testing that the output of the proposed Python tool works in this regard.

That would be really great! I'm a little worried about the NetCDF output from xarray, given @bradyrx's experience.

[mgduda]

I was able to run an idealized supercell simulation with MPAS-Atmosphere using a mesh generated by this script. For what it's worth, I installed xarray using pip.

[xylar]

@matthewhoffman, it sounds like you're in a minority in wanting the tool rewritten without xarray. Given that and my reluctance to do so, are you willing to approve the PR as it is or would you like to see further changes?

[matthewhoffman]

@xylar , if no one else objects to xarray, then I withdraw my concern. Feel free to merge this!

[xylar]

Thanks, @matthewhoffman. I appreciate that your goal was to keep the tools in this repo as clean an simple to use as possible. I'll definitely keep that goal in mind in my future additions.