Wrap-around the longitudes in the polygon locator #108
Conversation
Codecov ReportBase: 77.58% // Head: 77.56% // Decreases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## develop #108 +/- ##
===========================================
- Coverage 77.58% 77.56% -0.02%
===========================================
Files 805 805
Lines 57589 57606 +17
===========================================
+ Hits 44680 44683 +3
- Misses 12909 12923 +14
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
|
I think it is OK for now. @pmaciel probably has a better idea, but causing a lot of refactoring. That could be done in a follow up as needed. |
|
Yes @twsearle, there have been recent developments in eckit with an improved polygon class that correctly handles polygons spawning greater than 360 degree longitude. This is a starting point to introduce polygons as list of coordinates (the current case) and polygons as list of indices to a (fixed) list of coordinates -- such as a mesh element -- and working on (lat,lon), (x,y) and (x,y,z). These don't exist yet, but when they do they could support the partitioner as well as other fundamental functionality, too, also in other software in our stack (this updated polygon came requested from our software dowstream). There's some distance to cover still, however, so for the moment a solution like this works as well :-) |
|
Hi @pmaciel that sounds great! Yes I am happy with this interim solution for now. It would be good if the final solution can also be tested against the unusual orca edge case - a grid with a latitude domain ~ 60-420 or so for example! |
|
That's a good idea to test against that unusual case: would it be possible to contribute a test for this change, so that we can test when the transition to the improved polygon? The ORCA grids are important to us as well! |
|
@pmaciel, I have had a go at adding a test case in atlas-orca pr6, however in the course of doing that, I found another edge case I am not sure how solve. You can have a look there for details, but I think it is due to the particularities of how the atlas-orca partitioner is setup when dealing with the ORCA2 grids. These grids have a region that is offset in latitude longitude space, making a non-convex partition region. I think this is causing the partition polygon creation to fail. |
|
@twsearle, I'm not familiar with the ORCA-specific partitioner (@wdeconinck is). I can say, maybe stating the obvious, that partitioning elements (as opposed to points) does commonly generate "degenerate" polygons (where one vertex shares more than two edges) and less commonly disjoint polygons. Both these have to be considered by the partitioner, that it can handle these cases. This is caused by not just one of the grids proper (source or target) but sometimes the specific relation between them. It is still of great value to write a test, in this case more than one (right?), as precise as possible to trigger the issue, even if it isn't fixed right away. Maybe comment or disable the test out, but for me there's as much value in diagnosing the issue as with solving it ;-) |
|
This change breaks a nemo-feedback test: test_orcajedi_nemo_io_field_writer.x [1] �[31mTest "test parallel serially distributed write field array views" failed with unhandled eckit::Exception: Assertion failed: mesh_ in NodeColumns, line 209 of /home/h01/frwd/cylc-run/atlas108/share/mo-bundle/atlas/src/atlas/functionspace/NodeColumns.cc @ �[0m |
|
A quick update on the issue David noticed. We have confirmed that this is not associated with this PR and is in fact down to an issue with the test that occurs when atlas is updated from |
|
Hi @DJDavies2, @twsearle, please could you check if removing this assertion in "NodeColumns" fixes things? atlas/src/atlas/functionspace/NodeColumns.cc Line 209 in a63cff3 I'll merge this branch and the removal of this assertion in develop. |
* release/0.32.0: (24 commits) Version 0.32.0 Add function SolidBodyRotation Serial distribution to use entire mesh on every partition (#114) Adapt StructuredMeshGenerator to allow creation of meshes where some mpi tasks have no elements or nodes Allow interpolation to partitions without points or elements Remove assertion checking for empty mesh partition in NodeColumns (see #108) Wrap-around the longitudes in the polygon locator (#108) Github Actions: Prevent homebrew from auto-updating Github Actions: Prevent homebrew from auto-updating Github Actions: Update macos-11 to macos-12 (latest available) ATLAS-373 Enable gaussian latitudes computation for very high resolution, improving memory requirement with orders of magnitude Cosmetic change Add Fortran constructors for ShiftedLonLat, ShiftedLon, ShiftedLat grids ATLAS-372: Add missing header file Fixup previous commit, fixing segmentation fault Support more functionspaces for interpolation with FiniteElement, UnstructuredBilinearLonLat, GridBoxMethod, (K)NearestNeighbour Rename StructuredMeshGenerator option 3d -> three_dimensional, but still backward compatible Fix interpolation method UnstructuredBilinearLonLat for failing triangular elements Add missing CellColumns::halo() function Fix Mesh::operator bool() to reflect empty mesh ...
* release/0.32.0: (25 commits) Version 0.32.0 Fix compilation with ATLAS_BITS_LOCAL=64 Add function SolidBodyRotation Serial distribution to use entire mesh on every partition (#114) Adapt StructuredMeshGenerator to allow creation of meshes where some mpi tasks have no elements or nodes Allow interpolation to partitions without points or elements Remove assertion checking for empty mesh partition in NodeColumns (see #108) Wrap-around the longitudes in the polygon locator (#108) Github Actions: Prevent homebrew from auto-updating Github Actions: Prevent homebrew from auto-updating Github Actions: Update macos-11 to macos-12 (latest available) ATLAS-373 Enable gaussian latitudes computation for very high resolution, improving memory requirement with orders of magnitude Cosmetic change Add Fortran constructors for ShiftedLonLat, ShiftedLon, ShiftedLat grids ATLAS-372: Add missing header file Fixup previous commit, fixing segmentation fault Support more functionspaces for interpolation with FiniteElement, UnstructuredBilinearLonLat, GridBoxMethod, (K)NearestNeighbour Rename StructuredMeshGenerator option 3d -> three_dimensional, but still backward compatible Fix interpolation method UnstructuredBilinearLonLat for failing triangular elements Add missing CellColumns::halo() function ...
Allow the polyon locator to locate points in a grid when the xy space doesn't lie between 0-360 longitude.
JEDI uses ioda to distribute observations onto processors. This is done in ioda in AtlasDistribution to match the observation distribution to a chosen atlas grid. This uses the PolygonLocator to locate the owning partition for each observation.
This doesn't work on
atlas-orcagrids, because the xy space for orca data does not start at x=0 and doesn't finish at x=360 (see the meshgenerator).This PR is a proposal for how to solve this. When a point is not located by the locator, we try again but shifted by 360 degrees.
I cannot forsee how this could go wrong, however I don't know all the grids intimately, so I am very interested in your thoughts! I have rerun the tests with this change and things seem fine.
I have experimented with a test for the locator that catches the issue in atlas-orca and added it to atlas-orca PR#6, but this issue doesn't come up for other grids so I can't really cover it within atlas proper.