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

[odlomax]

Following on from PR #163, this PR primarily adds and tests the adjoint methods for the interpolation::method::SphericalVector class.

Several refactors were also applied to reduce complexity and code duplication:

• The matrix multiplication code was removed from SphericalVector.cc and given it's own class ComplexMatrixMultiply
• A SparseMatrix class was added that wraps the Eigen Sparse matrix into something that looks more like eckit::linalg::SparseMatrix with a Value template. This also confines all the nasty macros to one place.
• Types.h was added, containing all the main type definitions required for the SphericalVector class.
• ArrayForEachDim was added to array::helpers::ArrayForEach. This allows the iteration dimensions to be set using an std::integer_sequence.

Further testing was added

• Test added for ArrayForEachDim in ArrayForEach tests.
• Adjoint methods added to all interpolation tests.
• StructuredColumns to CubedSphere interpolation test added. Figures showing performance are given below. Note: StructuredColumns FixupHaloForVectors methods interferes with SphericalVector. You must make sure a halo exchange isn't called on the source field for it to work.

Notes:
Apologies for the strange commit history, I was trying to debug the MacOS build without a Mac.
Closes #166

Figures

O48 Structured Columns to Cubed Sphere interpolation error when using SphericalVector to treat vector fields (same setup as PR #163).

O48 Structured Columns to Cubed Sphere interpolation error when using FixupHaloForVectors to treat vector fields.

[odlomax]

@wdeconinck @pmaciel

Just thought you might be interested that I tried the vector interpolation scheme with StructuredColumns (see above). Also did some tests with double and half the source functionspace resolution. Reassuringly, there's quadratic convergence on the analytic solution. 😀

[wdeconinck]

Nice Addition @odlomax !

Could we try to split the PR in two, one for the arrayForEach and one for the adjoint?

We can merge the arrayForEach one first, and rebase this one after.

[odlomax]

Nice Addition @odlomax !

Could we try to split the PR in two, one for the arrayForEach and one for the adjoint?

We can merge the arrayForEach one first, and rebase this one after.

Absolutely. Watch this space...

[odlomax]

Nice Addition @odlomax !
Could we try to split the PR in two, one for the arrayForEach and one for the adjoint?
We can merge the arrayForEach one first, and rebase this one after.

Absolutely. Watch this space...

Done. Here's the PR.

[wdeconinck]

Other PR #171 is merged. Please rebase :)

[wdeconinck]

🥳 @odlomax

[odlomax]

🥳 @odlomax

Great success! 🤩

[odlomax]

No more pushes from me for now. I'll tidy up the macros when it's reviewed. 🙂

[wdeconinck]

How come you could drop this?

[odlomax]

I think the problem arose when supplying a template parameter functor to loop body. I've since separated out the use-cases into their own loops, as @pmaciel suggested. No longer triggers the error. 🙂

[wdeconinck]

I now played a bit with nvhpc/22.11 myself now in the routine ComplexMatrixMultiply::applyThreeVector
I tried simplifying until it worked, without using the arrayForeachDim construct. Following works:

atlas_omp_parallel_for( ... ) {
    ...
    if constexpr (Rank == 3) {
        for( idx_t jlev=0; jlev<targetSlice.shape(0); ++jlev) {
            const auto targetVector =
                complexWeight * Complex(sourceSlice(jlev,0), sourceSlice(jlev,1));

            targetSlice(jlev,0) += targetVector.real();
            targetSlice(jlev,1) += targetVector.imag();
            targetSlice(jlev,2) += realWeight * sourceSlice(jlev,2);
        }
    }
    if constexpr (Rank == 2) {
        ...
    }

What however doesn't work is creating a lambda that captures realWeight,complexWeight to do the inner part as:

atlas_omp_parallel_for( ... ) {
    ...
    auto matmul = [realWeight,complexWeight](auto&& sourceElem, auto&& targetElem) {
              const auto targetVector =
                  complexWeight * Complex(sourceElem(0), sourceElem(1));
              targetElem(0) += targetVector.real();
              targetElem(1) += targetVector.imag();
              targetElem(2) += realWeight * sourceElem(2);
    };

    if constexpr (Rank == 3) {
        for( idx_t jlev=0; jlev<targetSlice.shape(0); ++jlev) {
              auto sourceElem = sourceSlice.slice(jlev,array::Range::all());
              auto targetElem = targetSlice.slice(jlev,array::Range::all());
              matmul(sourceElem,targetElem);
        }
    }
    if constexpr (Rank == 2) {
        matmul(sourceSlice,targetSlice);
    }

If instead of "capturing" realWeight,complexWeight in the lambda, and passing it by argument, it WORKS again.

atlas_omp_parallel_for( ... ) {
    ...
    auto matmul = [](auto&& sourceElem, auto&& targetElem, auto realWeight, auto complexWeight) {
              const auto targetVector =
                  complexWeight * Complex(sourceElem(0), sourceElem(1));
              targetElem(0) += targetVector.real();
              targetElem(1) += targetVector.imag();
              targetElem(2) += realWeight * sourceElem(2);
    };

    if constexpr (Rank == 3) {
        for( idx_t jlev=0; jlev<targetSlice.shape(0); ++jlev) {
              auto sourceElem = sourceSlice.slice(jlev,array::Range::all());
              auto targetElem = targetSlice.slice(jlev,array::Range::all());
              matmul(sourceElem, targetElem, realWeight, complexWeight);
        }
    }
    if constexpr (Rank == 2) {
        matmul(sourceSlice, targetSlice, realWeight, complexWeight);
    }

So it appears the issue is with the lambda capture in the OpenMP region.
I verified that the problem is also still there with nvhpc/23.7

[odlomax]

I need to give up on templates. They're nothing but trouble! 😅

[wdeconinck]

I need to give up on templates. They're nothing but trouble! 😅

Sometimes it's really not worth the effort (probably many many many hours) to cater for added generality or save some lines of code. With if constexpr a lot of templates can be made a lot cleaner (and SFINAE can disappear).

[odlomax]

I need to give up on templates. They're nothing but trouble! 😅

Sometimes it's really not worth the effort (probably many many many hours) to cater for added generality or save some lines of code. With if constexpr a lot of templates can be made a lot cleaner (and SFINAE can disappear).

I'm glad concepts in C++20 has killed SFINAE for good. I think my problem is that lambda expressions have become my golden hammer, and now every problem looks like a nail!

[wdeconinck]

Thanks @odlomax really it's not your fault but nvidia compiler's fault. Thanks for working around! ❤️

[odlomax]

Thanks @odlomax really it's not your fault but nvidia compiler's fault. Thanks for working around! ❤️

Honestly, it's good life experience! We're going to have nvhpc on our new system, so it's useful to know its peculiarities!

[wdeconinck]

Looks good to me now. Do you mind if I "Squash-and-Merge" ?

[odlomax]

Looks good to me now. Do you mind if I "Squash-and-Merge" ?

Wonderful! Squash away!