Parallel transport generic vector interpolation method.

[odlomax]

Hi @wdeconinck,

I'm experimenting with a new method I've cooked up to perform geometry-agnostic vector interpolation. Essentially, it works out the angular difference between North and a great-circle arc connecting a source and a target point, then rotates the source field vector to match the target field vector.

The upshot of this is that we require four interpolation matrices instead of the usual one complex interpolation weights, but the results look promising. I've experimented below, interpolating from a CS-LFR-48 to an O48 (the interpolation method knows nothing about the cubed-sphere).

Here are some figures showing the results so far:

Fig 1: Eastward component of test field on CS-LFR-48 grid

Fig 2: Northward component of test field on CS-LFR-48 grid

Fig 3: Interpolation error when mapped on to O48 gird. Vector field components treated as independent scalars.

Fig 4: Interpolation error when mapped on to O48 gird. Source vectors rotated to match target vectors.

tl;dr I've made the error of interpolating winds over the poles go away, and I believe it will work for all matrix-based interpolation methods.

I'll try and get it finished up over the next week or so.

[codecov-commenter]

Codecov Report

Attention: 21 lines in your changes are missing coverage. Please review.

Comparison is base (24d175a) 79.92% compared to head (56badc7) 80.02%.

Files	Patch %	Lines
...polation/method/sphericalvector/SphericalVector.cc	86.71%	17 Missing :warning:
...rpolation/method/sphericalvector/SphericalVector.h	33.33%	2 Missing :warning:
src/atlas/util/Geometry.h	60.00%	2 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## develop #163 +/- ## =========================================== + Coverage 79.92% 80.02% +0.10% =========================================== Files 853 857 +4 Lines 63199 63498 +299 =========================================== + Hits 50511 50817 +306 + Misses 12688 12681 -7 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[wdeconinck]

Your updates to exclusively rely on Eigen is potentially problematic. Two ways...

1. Revert the change, which ensures that it will always work.
2. Check if compilation without Eigen works now (cmake : -DENABLE_EIGEN=OFF), and add some #if ATLAS_HAVE_EIGEN in SphericalVector.h and SphericalVector.cc to hide Eigen includes and types and then throw exceptions in the SphericalVector.cc when trying to use it. Also the compilation of SphericalVector should not be optional. The #if ATLAS_HAVE_EIGEN in code should take care of it.

[odlomax]

Your updates to exclusively rely on Eigen is potentially problematic. Two ways...

1. Revert the change, which ensures that it will always work.
2. Check if compilation without Eigen works now (cmake : -DENABLE_EIGEN=OFF), and add some #if ATLAS_HAVE_EIGEN in SphericalVector.h and SphericalVector.cc to hide Eigen includes and types and then throw exceptions in the SphericalVector.cc when trying to use it. Also the compilation of SphericalVector should not be optional. The #if ATLAS_HAVE_EIGEN in code should take care of it.

Thanks, @wdeconinck, I've replaced the optional compilation with pre-processor macros in relevant files. It builds and runs properly on my system with ENABLE_EIGEN set to both ON and OFF.

I did some digging in the current eckit sparse matrix. It does indeed prune out zero-valued triples (see here). This means all the real, imaginary and magnitude matrices can have different sparsity. The simplest in-code solution was to use Eigen directly. Minimal code changes would be required if the current eckit SparseMatrix is generalised to SparseMatrixT.

I'm going to open an issue which details what I think are the to-do items (hopefully) for follow-up PRs.

[wdeconinck]

I did some digging in the current eckit sparse matrix. It does indeed prune out zero-valued triples

That does seem like a sensible optimisation, but prevents us from fusing the two loops indeed, and then Eigen matrix solution is appropriate.

[odlomax]

I did some digging in the current eckit sparse matrix. It does indeed prune out zero-valued triples

That does seem like a sensible optimisation, but prevents us from fusing the two loops indeed, and then Eigen matrix solution is appropriate.

Its probably worth pointing out, I could fuse the the horizontal and vertical component loops, at the expense of adding some more complicated code.

[odlomax]

Hi @wdeconinck

I uncovered a lovely little feature where eckit and Eigen3 have slightly different sparsities in their CRS format. I've accounted for that now in the fused matrix-multiply loops. I've also ticked off the code-coverage tests in issue #166

I believe it's ready for review.

Cheers!

[pmaciel]

Hi @odlomax , what is the difference you found between Eigen's and eckit's CSR formats? I'm aware of a few (eckit doesn't require sorted column indices, not even uniqueness (which is a problem)) but is there any other difference? Your input might be invaluable here for an upcoming (separate) rewrite of the SparseMatrix.

I could only see a change re. the iterating through the compressed structure.

[odlomax]

Hi @odlomax , what is the difference you found between Eigen's and eckit's CSR formats? I'm aware of a few (eckit doesn't require sorted column indices, not even uniqueness (which is a problem)) but is there any other difference? Your input might be invaluable here for an upcoming (separate) rewrite of the SparseMatrix.

I could only see a change re. the iterating through the compressed structure.

On further consideration, my statement is not quite correct. The main difference is that Eigen's SparseMatrix constructor sorts each compressed row by j, whereas eckit does not.

The underlying data structures are the same, though. You can safely use the Eigen's Map template on an eckit SparseMatrix.

[wdeconinck]

intel/2021.4 failing on our HPC with an internal compiler error :(

 /usr/local/apps/intel/2021.4.0/compiler/latest/linux/bin/intel64/icpc -Datlas_EXPORTS -I../src -Isrc -I../atlas_io/src -Iatlas_io/src -I/usr/local/apps/fftw/3.3.10/INTEL/2021.4/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/fckit/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/fckit/module/fckit -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/fckit/include/fckit -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/geometry -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/linalg -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/maths -I/opt/mpi/openmpi/4.1.1.1/include -I/opt/mpi/openmpi/4.1.1.1/fortran-intel/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/mpi -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/option -I/usr/local/apps/eigen/3.4.0/include/eigen3 -O2 -g -DNDEBUG -fPIC -qopenmp -std=gnu++17 -MD -MT src/atlas/CMakeFiles/atlas.dir/interpolation/method/sphericalvector/SphericalVector.cc.o -MF src/atlas/CMakeFiles/atlas.dir/interpolation/method/sphericalvector/SphericalVector.cc.o.d -o src/atlas/CMakeFiles/atlas.dir/interpolation/method/sphericalvector/SphericalVector.cc.o -c ../src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
  ../src/atlas/interpolation/method/sphericalvector/SphericalVector.cc(70): internal error: assertion failed: find_assoc_pragma: pragma not found (il.c, line 25413 in find_assoc_pragma)

          functor(it.row(), it.col(), it.value());
          ^

  compilation aborted for ../src/atlas/interpolation/method/sphericalvector/SphericalVector.cc (code 4)

[odlomax]

intel/2021.4 failing on our HPC with an internal compiler error :(

 /usr/local/apps/intel/2021.4.0/compiler/latest/linux/bin/intel64/icpc -Datlas_EXPORTS -I../src -Isrc -I../atlas_io/src -Iatlas_io/src -I/usr/local/apps/fftw/3.3.10/INTEL/2021.4/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/fckit/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/fckit/module/fckit -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/fckit/include/fckit -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/geometry -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/linalg -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/maths -I/opt/mpi/openmpi/4.1.1.1/include -I/opt/mpi/openmpi/4.1.1.1/fortran-intel/include -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/mpi -I/etc/ecmwf/ssd/ssd1/tmpdirs/***.33879484/install/eckit/include/eckit/option -I/usr/local/apps/eigen/3.4.0/include/eigen3 -O2 -g -DNDEBUG -fPIC -qopenmp -std=gnu++17 -MD -MT src/atlas/CMakeFiles/atlas.dir/interpolation/method/sphericalvector/SphericalVector.cc.o -MF src/atlas/CMakeFiles/atlas.dir/interpolation/method/sphericalvector/SphericalVector.cc.o.d -o src/atlas/CMakeFiles/atlas.dir/interpolation/method/sphericalvector/SphericalVector.cc.o -c ../src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
  ../src/atlas/interpolation/method/sphericalvector/SphericalVector.cc(70): internal error: assertion failed: find_assoc_pragma: pragma not found (il.c, line 25413 in find_assoc_pragma)

          functor(it.row(), it.col(), it.value());
          ^

  compilation aborted for ../src/atlas/interpolation/method/sphericalvector/SphericalVector.cc (code 4)

Noooooo! 😢

A bit of Googling seems to suggest that the Intel compiler might be a bit fussy (buggy) when it comes to omp pragmas. It's also probable I've been a bit sloppy with my index types. If the latest push doesn't solve it, I'll rewrite the loops to operate over the SparseMatrix's internal C-arrays.

[wdeconinck]

Nice cleanup. Unfortunately still the same compile error with intel/2021.4 We could also disable OpenMP for that compiler and version as that thread you suggested says it could be fixed with 2022.2

[odlomax]

Nice cleanup. Unfortunately still the same compile error with intel/2021.4

Yeah, I just saw. What's confusing is that all the GitHub Actions pass on the JCSDA-internal fork end.

I guess I've got two options:

• Rewrite the function differently and hope for the best.
• Remove the OpenMP pragmas and accept that it's going to be slow for now (the error message is complaining about the pragma on the line above).

What would you suggest?

[odlomax]

Nice cleanup. Unfortunately still the same compile error with intel/2021.4 We could also disable OpenMP for that compiler and version as that thread you suggested says it could be fixed with 2022.2

It's certainly worth a try. Should at least help diagnose the problem! 😊

[odlomax]

Hi @odlomax,

as discussed this is a great contribution! As you see, I've only commented on technical details, because it looks like everything is well thought of. I admit, I think there is a bit overuse of templates to pass functors (when we really only handle the real and complex cases), but it is certainly not an important concern at the moment.

I've already put the "course" method as part of a revised eckit geometry library. I look forward to have ready the more flexible sparse matrix container and the sparse linear algebra backend too, as we discussed. It will happen in the coming months, and will also simplify much of the workarounds you had to painfuly go through.

I'm already aproving the PR -- this is ready to merge!

Thanks, @pmaciel, @wdeconinck,

I had a lot of fun with this one.

I freely admit I'm a template addict. If only we had C++ 20, I could at least throw a few concepts on them! I need to follow this up with an implementation of the adjoint methods. Shall I push @pmaciel's suggested changes to this PR, or the next one?

[pmaciel]

Thanks, @pmaciel, @wdeconinck,

I had a lot of fun with this one.

I freely admit I'm a template addict. If only we had C++ 20, I could at least throw a few concepts on them! I need to follow this up with an implementation of the adjoint methods. Shall I push @pmaciel's suggested changes to this PR, or the next one?

Since they're only technical would prefer that it's merged with the suggested changes :-)

[odlomax]

Thanks, @pmaciel, @wdeconinck, I had a lot of fun with this one. I freely admit I'm a template addict. If only we had C++ 20, I could at least throw a few concepts on them! I need to follow this up with an implementation of the adjoint methods. Shall I push @pmaciel's suggested changes to this PR, or the next one?

Since they're only technical would prefer that it's merged with the suggested changes :-)

It is done! Thanks, @pmaciel 🙂

[wdeconinck]

Thanks @odlomax for addressing further review. I will merge today when all tests have passed. This will be super useful.

[odlomax]

Thanks @odlomax for addressing further review. I will merge today when all tests have passed. This will be super useful.

Again, thanks @wdeconinck and @pmaciel. Always a pleasure.