Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add extended cubed-sphere dual halo with a correct alpha-beta projection. #90

Conversation

odlomax
Copy link
Contributor

@odlomax odlomax commented Feb 2, 2022

This is PRs 86 and 87 from JCSCA-internal, merged into one. The functionality added is as follows:

PR 86

This PR add functionality to the dual mesh halo. Now, the halo size can be set to any value < N-1, where N is the grid-size of a cubed sphere tile. This allows arbitrary-sized stencil operations to work on the cubed-sphere dual mesh.

An example of a cubed sphere dual mesh partition (8 PEs, default partitioner) with a halo of 3 is shown below.

dual-halo-2d
xy projection of an N48 partition with a halo of 3.

dual-halo-3d
xyz projection of an N48 partition with a halo of 3.

I've updated several of the existing cubed sphere dual mesh tests to use larger halo sizes.

PR 87

This PR predominantly adds an (x, y) to (alpha, beta) projection for tile t. Here, (alpha, beta) is defined strictly as the angular coordinates described by Ronchi et al., (1996).

In practice, the transform is a simple shift and rotate for xy coordinates internal to tile t. A non-linear stretch is applied for external points (i.e. halo xy coordinates).

The following example considers tile t=0 of the LFRic cubed sphere dual mesh, with halo = 3. The first image shows the xy projection. The second image gives the (alpha, beta) projection.

dual_halo_before
dual_halo_after

This required changes to several classes, summarised as follows:

  • grid::Tiles added tileCenter and tileJacobian which return the xy centre and d_xy/d_alphabeta Jacobian for each tile.
  • meshgenerator::detail::cubedsphere::Jacobian2 removed class and added equivalent to util namespace.
  • meshgenerator::detail::cubedsphere::NeighbourJacobian updated this class to use the new Tiles methods. This minor refactoring therefore serves as a test for the changes.
  • CubedSphereGrid added method to return a reference to the stored CubedSphereProjectionBase
  • redistribution::detail::RedistributionGeneric removed rank <= 3 constraint. Not related to projection, but it was an easy fix.

This transform is currently only implemented for the equiangular cubedsphere projection. It is also likely that I've introduced some unnecessary redundancy to the CubedSphereProjectionBase and Tiles classes. However, given that we intend to refactor the cubed sphere grid at some point, we'll eventually need to refactor these classes too.

@codecov-commenter
Copy link

codecov-commenter commented Feb 2, 2022

Codecov Report

Merging #90 (5afcfc8) into develop (290756e) will decrease coverage by 0.54%.
The diff coverage is 95.06%.

Impacted file tree graph

@@             Coverage Diff             @@
##           develop      #90      +/-   ##
===========================================
- Coverage    77.62%   77.08%   -0.55%     
===========================================
  Files          756      749       -7     
  Lines        51176    44935    -6241     
===========================================
- Hits         39727    34637    -5090     
+ Misses       11449    10298    -1151     
Impacted Files Coverage Δ
src/atlas/functionspace/CubedSphereColumns.h 100.00% <ø> (ø)
.../atlas/functionspace/detail/CubedSphereStructure.h 100.00% <ø> (ø)
src/atlas/grid/Tiles.h 100.00% <ø> (ø)
src/atlas/grid/detail/tiles/FV3Tiles.cc 89.62% <0.00%> (-5.75%) ⬇️
src/atlas/grid/detail/tiles/FV3Tiles.h 100.00% <ø> (ø)
src/atlas/grid/detail/tiles/LFRicTiles.h 100.00% <ø> (ø)
src/atlas/grid/detail/tiles/Tiles.cc 60.00% <ø> (-2.50%) ⬇️
...eshgenerator/detail/CubedSphereDualMeshGenerator.h 0.00% <ø> (-50.00%) ⬇️
...hgenerator/detail/cubedsphere/CubedSphereUtility.h 88.88% <ø> (-11.12%) ⬇️
.../projection/detail/CubedSphereEquiAnglProjection.h 71.42% <ø> (ø)
... and 541 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 290756e...5afcfc8. Read the comment docs.

@wdeconinck
Copy link
Member

Thanks @odlomax I will have look very soon

@odlomax
Copy link
Contributor Author

odlomax commented Feb 9, 2022

@wdeconinck, after our chat the other day, I removed JacobianXY class and added a general matrix class in atlas/util/Matrix.h. As per your suggestion, it uses eckit::maths::Matrix as a back end.

I also learned (the hard way) that composition > inheritance.

@odlomax
Copy link
Contributor Author

odlomax commented Feb 9, 2022

Hi @wdeconinck ,

I can't make sense of these linker errors when I try and build without Eigen:

/usr/bin/ld: lib/libatlas.so.0.27: undefined reference to `eckit::maths::lapack::getrf(int*, int*, double*, int*, int*, int*)'
  /usr/bin/ld: lib/libatlas.so.0.27: undefined reference to `eckit::maths::lapack::getri(int*, double*, int*, int*, double*, int*, int*)

Is it anything you've come across before?

@odlomax
Copy link
Contributor Author

odlomax commented Feb 9, 2022

It's also worth mentioning that transpose() in eckit/maths/MatrixLapack.h is non const when I think it should be const.

@odlomax
Copy link
Contributor Author

odlomax commented Feb 10, 2022

Hi @wdeconinck ,

I can't make sense of these linker errors when I try and build without Eigen:

/usr/bin/ld: lib/libatlas.so.0.27: undefined reference to `eckit::maths::lapack::getrf(int*, int*, double*, int*, int*, int*)'
  /usr/bin/ld: lib/libatlas.so.0.27: undefined reference to `eckit::maths::lapack::getri(int*, double*, int*, int*, double*, int*, int*)

Is it anything you've come across before?

My bad. Fixed it!

@odlomax odlomax closed this Feb 10, 2022
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