Parallel transport generic vector interpolation method. #163

odlomax · 2023-11-17T17:02:41Z

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.

make_view calls erroneously in loop body.

codecov-commenter · 2023-11-30T09:18:24Z

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 ⚠️
...rpolation/method/sphericalvector/SphericalVector.h	33.33%	2 Missing ⚠️
src/atlas/util/Geometry.h	60.00%	2 Missing ⚠️

Additional details and impacted files

@@             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

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

wdeconinck · 2023-11-30T09:32:38Z

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

Revert the change, which ensures that it will always work.
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.

…e and header.

odlomax · 2023-11-30T11:56:26Z

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

Revert the change, which ensures that it will always work.

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.

src/atlas/interpolation/method/MethodFactory.cc

wdeconinck · 2023-11-30T13:32:10Z

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 · 2023-11-30T13:44:39Z

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.

wdeconinck · 2023-11-30T15:22:11Z

src/atlas/interpolation/method/MethodFactory.cc

It was still good to have it in this file, but not guarded by #if ATLAS_HAVE_EIGEN if that makes sense.

I can see that it adds unnecessary marcro mess. I'm just marvelling at how the factories can register at run time instead of compile time!

What I meant is what I added to your branch with a96057e

Ah! That looks far safer!

… Write test for 3-vector field (note: this works when hacking current 2-vector test).

added tests for 3-vector and 2 vector fields.

wdeconinck · 2023-12-08T16:42:52Z

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 · 2023-12-08T16:54:21Z

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 · 2023-12-08T17:00:57Z

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 · 2023-12-12T13:54:02Z

src/atlas/interpolation/method/sphericalvector/SphericalVector.cc

+#define atlas_omp_parallel_for for
+#endif
+#endif
+


Thanks, @wdeconinck. Hit it with macros until the problem goes away! ❤

The problem indeed went away 😃

wdeconinck · 2023-12-13T16:10:08Z

src/atlas/interpolation/method/sphericalvector/SphericalVector.cc

+    const auto sourceVector = Complex(sourceVars(0), sourceVars(1));
+    const auto targetVector = complexWeight * sourceVector;
+    targetVars(0) += targetVector.real();
+    targetVars(1) += targetVector.imag();


Without going into the detailed maths, why do you have to
targetView.assign(0)
and then use += for the targetVars ? Can this not be replaced with just = ?

Never mind, I understand, due to the matrix multiply internally.
But perhaps performance can be gained if the assign(0) could be removed and done inside the matrixMultiply

It's because each target element (named targetVars here) has a contribution from multiple source elements (sourceVars). If I changed += to =, it would set the target element to the last source element on that row of the matrix.

I could do the assign(0.) within matrixMultiply function, but it would have to be on the first line, and therefore be the equivalent operation as it currently is. The SparseMatrixForEach is the only function that has knowledge of i as a whole row, but it doesn't know what functor is doing (it's used in a couple of different contexts).

I suspect details like this will have to be refactored to remove code duplication when the adjoint methods are added. In that context, you don't assign zero to the array you're writing to.

OK it can stay as-is. We should not attempt premature optimisation now.

pmaciel

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!

src/atlas/interpolation/method/sphericalvector/SphericalVector.h

src/atlas/interpolation/method/sphericalvector/SphericalVector.cc

pmaciel · 2023-12-14T16:18:13Z

src/atlas/interpolation/method/sphericalvector/SphericalVector.cc

+    return;
+  }
+
+  ATLAS_ASSERT_MSG(sourceField.variables() == 2 || sourceField.variables() == 3,


Shouldn't this assertion come right at the beggining of the method?

I'd also have to check that the field was a vector type. By this point in the function, it definitely is.

odlomax · 2023-12-14T17:03:21Z

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 · 2023-12-14T17:15:41Z

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 :-)

…-internal/atlas into feature/parallel_transport

odlomax · 2023-12-14T18:39:52Z

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 · 2023-12-15T07:24:15Z

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

odlomax · 2023-12-15T11:46:20Z

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.

Implemented forward interpolation. TODO: adjoint and more tests.

2c9512f

github-actions bot added the contributor label Nov 17, 2023

odlomax added 9 commits November 20, 2023 18:36

Fixed race condition errors.

7c0a5d7

Replaced rotation matrices with complex interpolation weights.

d01f5dd

Merge branch 'develop' into feature/parallel_transport

61ec2e7

Fixed typos before discussion.

d5c9ef5

Renamed class to SphericalVector. Added Eigen3 matrices.

088403e

Update SphericalVector.cc

a6d0df9

make_view calls erroneously in loop body.

Tidied up SphericalVector class.

5064f44

Finalising for PR.

d0da681

Merge branch 'develop' into feature/parallel_transport

91de6bf

wdeconinck added the approved-for-ci label Nov 30, 2023

Minor cosmetic changes.

6b43f84

github-actions bot removed the approved-for-ci label Nov 30, 2023

odlomax added 2 commits November 30, 2023 11:33

Replaced optional compilation with #if ATLAS_HAVE_EIGEN in source fil…

35611a8

…e and header.

Removed redundant macros.

a913b18

odlomax marked this pull request as ready for review November 30, 2023 11:48

odlomax mentioned this pull request Nov 30, 2023

Adding robust vector-field interpolation functionality (current PR and further work). #166

Closed

6 tasks

wdeconinck reviewed Nov 30, 2023

View reviewed changes

src/atlas/interpolation/method/MethodFactory.cc Outdated Show resolved Hide resolved

Removed static factory linking.

c4481aa

wdeconinck reviewed Nov 30, 2023

View reviewed changes

odlomax added 3 commits December 1, 2023 18:28

Fused horizontal and vertical component matrix-multiplications. TODO:…

6107e91

… Write test for 3-vector field (note: this works when hacking current 2-vector test).

Tidied fused loop.

77856f3

Uncovered and fixed differences in eckit and Eigen3 CRS format. Also

346f7f1

added tests for 3-vector and 2 vector fields.

wdeconinck added the approved-for-ci label Dec 5, 2023

wdeconinck added the approved-for-ci label Dec 8, 2023

Merge branch 'develop' into feature/parallel_transport

ba8193c

github-actions bot removed the approved-for-ci label Dec 11, 2023

Disable OpenMP for older intel-classic compiler (< intel/2022.2)

fd77d76

wdeconinck added the approved-for-ci label Dec 12, 2023

odlomax commented Dec 12, 2023

View reviewed changes

wdeconinck added 2 commits December 12, 2023 15:25

Enable test with CONDITION statement

3686d84

Revert whitespace changes

09af7e7

github-actions bot removed the approved-for-ci label Dec 12, 2023

wdeconinck added 2 commits December 12, 2023 16:39

Make greatCircleCourse private before moving to eckit

1212789

Fix header includes

98ab54a

wdeconinck reviewed Dec 13, 2023

View reviewed changes

wdeconinck approved these changes Dec 13, 2023

View reviewed changes

wdeconinck added the approved-for-ci label Dec 13, 2023

pmaciel approved these changes Dec 14, 2023

View reviewed changes

odlomax added 2 commits December 14, 2023 17:44

Addressed reviewer comments.

96b28b2

Merge branch 'feature/parallel_transport' of https://github.com/JCSDA…

56badc7

…-internal/atlas into feature/parallel_transport

github-actions bot removed the approved-for-ci label Dec 14, 2023

wdeconinck added the approved-for-ci label Dec 15, 2023

wdeconinck merged commit e734e7a into ecmwf:develop Dec 15, 2023
92 checks passed

odlomax deleted the feature/parallel_transport branch January 2, 2024 11:04

odlomax mentioned this pull request Jan 8, 2024

Refactoring of interpolation::method::SphericalVector and implementation of adjoint methods. #168

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Parallel transport generic vector interpolation method. #163

Parallel transport generic vector interpolation method. #163

odlomax commented Nov 17, 2023 •

edited

codecov-commenter commented Nov 30, 2023 •

edited

wdeconinck commented Nov 30, 2023 •

edited

odlomax commented Nov 30, 2023

wdeconinck commented Nov 30, 2023

odlomax commented Nov 30, 2023

wdeconinck Nov 30, 2023

odlomax Nov 30, 2023 •

edited

wdeconinck Dec 7, 2023

odlomax Dec 7, 2023

wdeconinck commented Dec 8, 2023 •

edited

odlomax commented Dec 8, 2023

odlomax commented Dec 8, 2023

odlomax Dec 12, 2023

wdeconinck Dec 12, 2023

wdeconinck Dec 13, 2023

wdeconinck Dec 13, 2023

odlomax Dec 13, 2023

odlomax Dec 13, 2023

wdeconinck Dec 13, 2023

pmaciel left a comment

pmaciel Dec 14, 2023

odlomax Dec 14, 2023

odlomax commented Dec 14, 2023

pmaciel commented Dec 14, 2023

odlomax commented Dec 14, 2023

wdeconinck commented Dec 15, 2023

odlomax commented Dec 15, 2023

Parallel transport generic vector interpolation method. #163

Parallel transport generic vector interpolation method. #163

Conversation

odlomax commented Nov 17, 2023 • edited

codecov-commenter commented Nov 30, 2023 • edited

Codecov Report

wdeconinck commented Nov 30, 2023 • edited

odlomax commented Nov 30, 2023

wdeconinck commented Nov 30, 2023

odlomax commented Nov 30, 2023

Choose a reason for hiding this comment

odlomax Nov 30, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wdeconinck commented Dec 8, 2023 • edited

odlomax commented Dec 8, 2023

odlomax commented Dec 8, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pmaciel left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

odlomax commented Dec 14, 2023

pmaciel commented Dec 14, 2023

odlomax commented Dec 14, 2023

wdeconinck commented Dec 15, 2023

odlomax commented Dec 15, 2023

odlomax commented Nov 17, 2023 •

edited

codecov-commenter commented Nov 30, 2023 •

edited

wdeconinck commented Nov 30, 2023 •

edited

odlomax Nov 30, 2023 •

edited

wdeconinck commented Dec 8, 2023 •

edited