Unstructured MPAS meshes cannot be used with conservative, bilinear, and patch mapping files in a compset

[douglasjacobsen]

Currently, when creating mapping files for an MPAS mesh that has cells removed (i.e. land for the ocean / sea ice) some cells are not map-able (according to ESMF) for each mapping type (conservative, bilinear, and patch).

Typically, we run a simulation where fluxes are remapped using conservative, state fields are remapped using bilinear, and vectors are remapped using patch.

When creating domain files from each of these different mapping files, the mask of cells to be mapped is different, meaning a single domain file cannot be used for all three of the mapping files. In this case, the bilinear / patch mapping files have less ocean cells that are mapped to than the conservative does.

This causes two issues:

1. When running with the conservative domain file, and the bilinear / patch mapping files the coupler sends garbage to the ocean when trying to map to a cell that has no remapping weights. These end up killing the model, since things like wind stress are remapped and given random values (some of which are of the order 1e30).
2. When running with the bilinear / patch domain file and the conservative mapping file, the coupler will throw away some fluxes, causing the resulting simulation to be less conservative than it should be.

I can think of four solutions, none of which I know how to carry out:

1. Within the coupler, if a cell should be mapped to, but doesn't have weights, remap with weights of 0.
2. Within ESMF, if a cell cannot be mapped, give it weights of 0
3. Within ESMF, if a cell cannot be mapped, give it weights equivalent to one of it's neighbors (that are non-zero).
4. Within ESMF, ensure all cells are given correct remapping weights

In addition to this, I have no idea how to test a mapping file. It would be good to know if there is something like taking an atmospheric field of 1s everywhere, and applying a remapping to it, and making sure the values on the ocean grid are all 1s (or some appropriate value).

[mt5555]

one point of clarification: When making maps for use in ACME, at least for the CAM-SE/POP case, all the maps are made without a mask, and are global grids with no points removed. So are you talking about this first step - in which case maybe you should not be removing cells from the MPAS mesh? Or are you talking about the second step - building the domain files, which takes as input the ESMF maps without the last mask?

[mt5555]

And for the record, I believe both state variables and fluxes have to be mapped with the same map file when going from ocean -> atmosphere. I think this constraint is related to this issue. This is the problematic direction because the land/ocean mask is defined on the ocean grid, and then the mapped version is imposed on the atmosphere grid.

[douglasjacobsen]

@mt5555 POP experience doesn't really help with MPAS. We don't have a mesh we run with the contains all the cells, so we can't create mapping files for an "unmasked" ocean grid.

In POP, all cells exist all the time, so the mesh they run with is the same as the mesh they create the mapping files with.

And, I'd argue for performance reasons, not removing the cells from the MPAS mesh is not a viable solution.

Currently (as I understand it at least, but I could be wrong), atmosphere -> ocean fluxes are conservative, but state fields use bilinear, and vectors use patch.

[mt5555]

So if the MPAS grid has "holes" in it, I think you will be in uncharted territory with respect to ACME/CESM and the coupler!

Agreed: atmosphere->ocean, the choice of mapping files is not constrained. But ocean->atmosphere, all mapping files have to be the same, due to how the land mask is implemented on the atmosphere grid. I think if you can fix the issue you describe above, the same approach may also allow us to remove this ocn->atm map file constraint.

[douglasjacobsen]

@mt5555 My main issue is that I don't know how to go about fixing it.

I'm not sure where to look in the coupler, and I definitely don't know anything about fixing issues within ESMF.

[mt5555]

How about, as an interim solution, making all mapping files conservative?

[douglasjacobsen]

@mt5555 That's what we're currently doing for the C compset runs. That works fine, but it gives pretty crappy results. So, this "bug" is more something that needs to be fixed before we can do semi-realistic simulations on the November-ish time frame.

[maltrud]

as a point of clarification: the reason this problem showed up in the first place is because the esmf mapping file generator fails to assign weights for 257 ocean cells for the EC60-30km mpas grid when doing bilinear and patch (all ocean cells get weights for conservative). there is an explanation for this, which Doug understands much better than i do. the reason we're thinking that the coupler may need to be modified is because we have no leverage with esmf to fix the problem. i personally don't like the idea of having to zero out some of the forcing fields at 257 locations do to problems with generating remapping weights.

[mt5555]

Have you tried the original SCRIP code- it's difficult to use for conservative maps, but its bilinear and distance weighted maps are quite robust.

Given a global atmospheric grid as the source grid, I would think both bilinear and patch maps should be able to map to every single point in the target grid. If Doug has the time, I would be interested in understanding why ESMF could not map to these 257 points. It seem like a bug in ESMF (hence the suggestion to try SCRIP)

[maltrud]

i suggested using, but scrip doesn't support it since bilinear isn't strictly defined going to a grid that doesn't have 4 sides. i think esmf does something using triangles which isn't exactly 'bilinear' either, but close enough i guess. i also suggested hacking something into scrip that we could use but that was met with a negative response since then it isn't consistent with acme principles of using one-off software.

i was thinking maybe we could ask esmf to use nearest neighbor in those cases where patch and bilinear fail. at least we'd have something realistic there.

[mt5555]

I interpreted above that the problem was atmosphere->ocean maps having points on the ocean grid that no data was mapped into them. So for bilinear, only the atmosphere grid needs to have 4 points containing every ocean point, which will always be the case?

Or is the problem the other direction, ocean->atmosphere? In which case SCRIP nearest neighbor should give weights to every ocean grid point.

[maltrud]

you're right, Mark. i had it backwards--there is no reason why esmf should have a problem going from atm=>ocn using bilinear since it is the source grid that needs to have 4 corners. so i took your advice and tried using scrip and it worked fine. i'm not sure why it failed the first time i tried this. most likely my brain failed. i have a job waiting to test the new remapping file.

interesting that there was no apparent trouble for esmf going from ocn=>atm. of course, when i ran check_maps on these, they reported passing on all 14 tests, though we now know that they weren't ok. but i guess in an RMS sense, 257 points won't affect the global metric.

assuming this works, we will face the issue of scrip not having a 'patch' option. i'm sure we (ie, Phil J, the grand poobah of scrip) could add it pretty easily.

[mt5555]

I've been testing map files generated with Paul Ullrich's TempestRemap. As part of this work, I've been comparing SE->FV maps generated by ESMF and TempestRemap. Some things of interested to this github issue:

(1) I've found the same problem with missing grid points in the ESMF bilinear and patch based maps, but have not yet seen any missing points with TempestRemap maps. So switching to TempestRemap maps may solve this problem.

(2) ESMF "patch", for a variety of analytic test functions, is no more accurate than bilinear.

(3) TempestRemap conservative non-monotone map is significantly more accurate than patch, bilinear or aave. It's probably the best map to use for mapping velocity fields

(4) TempestRemap conserve+mono is slightly more accurate that 'aave' maps.

[douglasjacobsen]

@mt5555: Thanks for the update. How does one use TempestRemap?

Aside from that ESMF has released a new beta snapshot (53 specifically) that fixes the missing point issue with maps from ATM->OCN, but not the other way around (at least when OCN is MPAS in this case).

[paullric]

Hi Doug,

The TempestRemap repository is here: https://github.com/ClimateGlobalChange/tempestremap

Basically it requires either an Exodus or SCRIP file describing the source and target meshes. An overlap mesh (the intersection of the source and target meshes) is then created using GenerateOverlapMesh and the mapping weights generated via GenerateOfflineMap. The command-line options are mostly documented in the README.md file.

Note that the use of an incomplete target mesh has been tested and seems to work well, but I haven't tried the weight generator with an incomplete source mesh. If you could provide me a low-resolution (25km or coarser) MPAS ocean mesh, I can validate it is working correctly and provide the command line sequence necessary for generating the maps.

[douglasjacobsen]

@paullric: You can find a 120km MPAS-Ocean SCRIP file here to play with: https://acme-svn2.ornl.gov/acme-repo/acme/inputdata/share/domains/domain.ocn.mpas120.121116.nc

Thanks for all the information, I'll let you know if I run into any issues.

[rljacob]

Can this be closed?

[douglasjacobsen]

I think the answer is no. As far as I know both tools still have issues creating mapping files when both meshes have holes in them (which we plan on doing for ocean <-> land ice maps).

That being said, if you'd prefer we could close this issue, and make a new issues related to those other mapping files, but a lot of the information for that other issue is contained in this one.