From 5aab9e2e75ce618651736b91648619ef9645f5d6 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 13 Feb 2023 11:34:09 +0100 Subject: [PATCH 01/19] C++17 is now a requirement --- CMakeLists.txt | 1 - atlas_io/CMakeLists.txt | 6 +----- atlas_io/src/atlas_io/CMakeLists.txt | 4 +--- cmake/atlas_compile_flags.cmake | 2 +- cmake/features/CXX17.cmake | 5 ----- 5 files changed, 3 insertions(+), 15 deletions(-) delete mode 100644 cmake/features/CXX17.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index bbb27b8a8..d44f51625 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,7 +55,6 @@ include( features/INCLUDE_WHAT_YOU_USE ) include( features/INIT_SNAN ) include( features/DOCS ) include( features/ATLAS_RUN ) -include( features/CXX17 ) ################################################################################ # sources diff --git a/atlas_io/CMakeLists.txt b/atlas_io/CMakeLists.txt index cf29b9a4d..26f549808 100644 --- a/atlas_io/CMakeLists.txt +++ b/atlas_io/CMakeLists.txt @@ -28,14 +28,10 @@ ecbuild_debug( " eckit_FEATURES : [${eckit_FEATURES}]" ) ################################################################################ # Features that can be enabled / disabled with -DENABLE_ -ecbuild_add_option( FEATURE CXX17 - DESCRIPTION "Use C++17 standard" - DEFAULT OFF ) - ################################################################################ # sources -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) check_cxx_source_compiles( "#include \n int main() { char * type; int status; char * r = abi::__cxa_demangle(type, 0, 0, &status); }" diff --git a/atlas_io/src/atlas_io/CMakeLists.txt b/atlas_io/src/atlas_io/CMakeLists.txt index 97140b00a..395d827d2 100644 --- a/atlas_io/src/atlas_io/CMakeLists.txt +++ b/atlas_io/src/atlas_io/CMakeLists.txt @@ -91,6 +91,4 @@ ecbuild_add_library( TARGET atlas_io ${CMAKE_CURRENT_BINARY_DIR}/detail/defines.h ) -target_compile_features( atlas_io PUBLIC - cxx_std_11 - $<${HAVE_CXX17}:cxx_std_17> ) +target_compile_features( atlas_io PUBLIC cxx_std_17 ) diff --git a/cmake/atlas_compile_flags.cmake b/cmake/atlas_compile_flags.cmake index 235f12668..9b350d478 100644 --- a/cmake/atlas_compile_flags.cmake +++ b/cmake/atlas_compile_flags.cmake @@ -6,7 +6,7 @@ # granted to it by virtue of its status as an intergovernmental organisation nor # does it submit to any jurisdiction. -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) if( CMAKE_CXX_COMPILER_ID MATCHES Cray ) diff --git a/cmake/features/CXX17.cmake b/cmake/features/CXX17.cmake deleted file mode 100644 index 7c6864031..000000000 --- a/cmake/features/CXX17.cmake +++ /dev/null @@ -1,5 +0,0 @@ -### C++17 ... - -ecbuild_add_option( FEATURE CXX17 - DESCRIPTION "Use C++17 standard" - DEFAULT OFF ) From 3f61cbb4d2a8e3540321745e3f548028298d4f78 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 14 Feb 2023 11:11:23 +0100 Subject: [PATCH 02/19] Modernise, using make_unique, and auto return type for detail return types --- src/atlas/array/gridtools/GridToolsArray.cc | 6 ++--- src/atlas/array/native/NativeArray.cc | 18 +++++++-------- src/atlas/array/native/NativeArrayView.h | 10 --------- src/atlas/grid/detail/grid/CubedSphere.h | 14 ++++++------ src/atlas/grid/detail/grid/Structured.h | 8 +++---- src/atlas/grid/detail/grid/Unstructured.h | 8 +++---- src/atlas/library/Library.cc | 2 +- src/atlas/util/Bitflags.h | 25 +++++++++------------ 8 files changed, 39 insertions(+), 52 deletions(-) diff --git a/src/atlas/array/gridtools/GridToolsArray.cc b/src/atlas/array/gridtools/GridToolsArray.cc index 27a5ee462..5c5a0c5cb 100644 --- a/src/atlas/array/gridtools/GridToolsArray.cc +++ b/src/atlas/array/gridtools/GridToolsArray.cc @@ -67,7 +67,7 @@ class ArrayT_impl { static_assert(sizeof...(UInts) > 0, "1"); auto gt_storage = create_gt_storage::type>(dims...); using data_store_t = typename std::remove_pointer::type; - array_.data_store_ = std::unique_ptr(new GridToolsDataStore(gt_storage)); + array_.data_store_ = std::unique_ptr>(gt_storage); array_.spec_ = make_spec(gt_storage, dims...); } @@ -77,7 +77,7 @@ class ArrayT_impl { auto gt_storage = create_gt_storage::type, ::gridtools::alignment>(dims...); using data_store_t = typename std::remove_pointer::type; - array_.data_store_ = std::unique_ptr(new GridToolsDataStore(gt_storage)); + array_.data_store_ = std::make_unique>(gt_storage); array_.spec_ = make_spec(gt_storage, dims...); } @@ -230,7 +230,7 @@ class ArrayT_impl { void construct_with_layout(UInts... dims) { auto gt_data_store_ptr = create_gt_storage(dims...); using data_store_t = typename std::remove_pointer::type; - array_.data_store_ = std::unique_ptr(new GridToolsDataStore(gt_data_store_ptr)); + array_.data_store_ = std::make_unique>(gt_data_store_ptr); array_.spec_ = make_spec(gt_data_store_ptr, dims...); } diff --git a/src/atlas/array/native/NativeArray.cc b/src/atlas/array/native/NativeArray.cc index e75defeae..015f1f521 100644 --- a/src/atlas/array/native/NativeArray.cc +++ b/src/atlas/array/native/NativeArray.cc @@ -116,27 +116,27 @@ ArrayT::ArrayT(ArrayDataStore* ds, const ArraySpec& spec) { template ArrayT::ArrayT(idx_t dim0) { spec_ = ArraySpec(make_shape(dim0)); - data_store_ = std::unique_ptr(new native::DataStore(spec_.size())); + data_store_ = std::make_unique>(spec_.size()); } template ArrayT::ArrayT(idx_t dim0, idx_t dim1) { spec_ = ArraySpec(make_shape(dim0, dim1)); - data_store_ = std::unique_ptr(new native::DataStore(spec_.size())); + data_store_ = std::make_unique>(spec_.size()); } template ArrayT::ArrayT(idx_t dim0, idx_t dim1, idx_t dim2) { spec_ = ArraySpec(make_shape(dim0, dim1, dim2)); - data_store_ = std::unique_ptr(new native::DataStore(spec_.size())); + data_store_ = std::make_unique>(spec_.size()); } template ArrayT::ArrayT(idx_t dim0, idx_t dim1, idx_t dim2, idx_t dim3) { spec_ = ArraySpec(make_shape(dim0, dim1, dim2, dim3)); - data_store_ = std::unique_ptr(new native::DataStore(spec_.size())); + data_store_ = std::make_unique>(spec_.size()); } template ArrayT::ArrayT(idx_t dim0, idx_t dim1, idx_t dim2, idx_t dim3, idx_t dim4) { spec_ = ArraySpec(make_shape(dim0, dim1, dim2, dim3, dim4)); - data_store_ = std::unique_ptr(new native::DataStore(spec_.size())); + data_store_ = std::make_unique>(spec_.size()); } template @@ -146,20 +146,20 @@ ArrayT::ArrayT(const ArrayShape& shape) { for (size_t j = 0; j < shape.size(); ++j) { size *= size_t(shape[j]); } - data_store_ = std::unique_ptr(new native::DataStore(size)); + data_store_ = std::make_unique>(size); spec_ = ArraySpec(shape); } template ArrayT::ArrayT(const ArrayShape& shape, const ArrayAlignment& alignment) { spec_ = ArraySpec(shape, alignment); - data_store_ = std::unique_ptr(new native::DataStore(spec_.allocatedSize())); + data_store_ = std::make_unique>(spec_.allocatedSize()); } template ArrayT::ArrayT(const ArrayShape& shape, const ArrayLayout& layout) { spec_ = ArraySpec(shape); - data_store_ = std::unique_ptr(new native::DataStore(spec_.size())); + data_store_ = std::make_unique>(spec_.size()); for (size_t j = 0; j < layout.size(); ++j) { ATLAS_ASSERT(spec_.layout()[j] == layout[j]); } @@ -167,7 +167,7 @@ ArrayT::ArrayT(const ArrayShape& shape, const ArrayLayout& layout) { template ArrayT::ArrayT(ArraySpec&& spec): Array(std::move(spec)) { - data_store_ = std::unique_ptr(new native::DataStore(spec_.allocatedSize())); + data_store_ = std::make_unique>(spec_.allocatedSize()); } template diff --git a/src/atlas/array/native/NativeArrayView.h b/src/atlas/array/native/NativeArrayView.h index 07c436024..b902cbb95 100644 --- a/src/atlas/array/native/NativeArrayView.h +++ b/src/atlas/array/native/NativeArrayView.h @@ -291,24 +291,14 @@ class ArrayView { /// auto slice3 = view.slice( Range::all(), Range::all(), Range::dummy() ); /// @endcode template -#ifndef DOXYGEN_SHOULD_SKIP_THIS - typename slice_t::type slice(Args... args) { -#else - // C++14 will allow auto return type auto slice(Args... args) { -#endif return slicer_t(*this).apply(args...); } /// @brief Obtain a slice from this view: view.slice( Range, Range, ... ) template -#ifndef DOXYGEN_SHOULD_SKIP_THIS - typename const_slice_t::type slice(Args... args) const { -#else - // C++14 will allow auto return type auto slice(Args... args) const { -#endif return const_slicer_t(*this).apply(args...); } diff --git a/src/atlas/grid/detail/grid/CubedSphere.h b/src/atlas/grid/detail/grid/CubedSphere.h index 998bc30be..69d29b004 100644 --- a/src/atlas/grid/detail/grid/CubedSphere.h +++ b/src/atlas/grid/detail/grid/CubedSphere.h @@ -345,7 +345,7 @@ class CubedSphere : public Grid { // Note that i is the fastest index, followed by j, followed by t std::unique_ptr nextElement(const idx_t i, const idx_t j, const idx_t t) const { - std::unique_ptr ijt(new int[3]); + auto ijt = std::make_unique(3); ijt[0] = i; ijt[1] = j; @@ -396,22 +396,22 @@ class CubedSphere : public Grid { // ---------------------------- virtual std::unique_ptr xy_begin() const override { - return std::unique_ptr(new IteratorXY(*this)); + return std::make_unique(*this); } virtual std::unique_ptr xy_end() const override { - return std::unique_ptr(new IteratorXY(*this, false)); + return std::make_unique(*this, false); } virtual std::unique_ptr lonlat_begin() const override { - return std::unique_ptr(new IteratorLonLat(*this)); + return std::make_unique(*this); } virtual std::unique_ptr lonlat_end() const override { - return std::unique_ptr(new IteratorLonLat(*this, false)); + return std::make_unique(*this, false); } virtual std::unique_ptr tij_begin() const { - return std::unique_ptr(new IteratorTIJ(*this)); + return std::make_unique(*this); } virtual std::unique_ptr tij_end() const { - return std::unique_ptr(new IteratorTIJ(*this, false)); + return std::make_unique(*this, false); } // Default configurations diff --git a/src/atlas/grid/detail/grid/Structured.h b/src/atlas/grid/detail/grid/Structured.h index 2783c022a..bd3687528 100644 --- a/src/atlas/grid/detail/grid/Structured.h +++ b/src/atlas/grid/detail/grid/Structured.h @@ -337,16 +337,16 @@ class Structured : public Grid { const YSpace& yspace() const { return yspace_; } virtual std::unique_ptr xy_begin() const override { - return std::unique_ptr(new IteratorXY(*this)); + return std::make_unique(*this); } virtual std::unique_ptr xy_end() const override { - return std::unique_ptr(new IteratorXY(*this, false)); + return std::make_unique(*this, false); } virtual std::unique_ptr lonlat_begin() const override { - return std::unique_ptr(new IteratorLonLat(*this)); + return std::make_unique(*this); } virtual std::unique_ptr lonlat_end() const override { - return std::unique_ptr(new IteratorLonLat(*this, false)); + return std::make_unique(*this, false); } gidx_t index(idx_t i, idx_t j) const { return jglooff_[j] + i; } diff --git a/src/atlas/grid/detail/grid/Unstructured.h b/src/atlas/grid/detail/grid/Unstructured.h index 30b09b995..e5a07848d 100644 --- a/src/atlas/grid/detail/grid/Unstructured.h +++ b/src/atlas/grid/detail/grid/Unstructured.h @@ -175,16 +175,16 @@ class Unstructured : public Grid { } virtual std::unique_ptr xy_begin() const override { - return std::unique_ptr(new IteratorXY(*this)); + return std::make_unique(*this); } virtual std::unique_ptr xy_end() const override { - return std::unique_ptr(new IteratorXY(*this, false)); + return std::make_unique(*this, false); } virtual std::unique_ptr lonlat_begin() const override { - return std::unique_ptr(new IteratorLonLat(*this)); + return std::make_unique(*this); } virtual std::unique_ptr lonlat_end() const override { - return std::unique_ptr(new IteratorLonLat(*this, false)); + return std::make_unique(*this, false); } Config meshgenerator() const override; diff --git a/src/atlas/library/Library.cc b/src/atlas/library/Library.cc index 10c06063d..0ce439bb6 100644 --- a/src/atlas/library/Library.cc +++ b/src/atlas/library/Library.cc @@ -305,7 +305,7 @@ void Library::initialise(const eckit::Parametrisation& config) { Adaptor(const eckit::CodeLocation& loc, const std::string& title): trace{loc, title} {} atlas::Trace trace; }; - return std::unique_ptr(new Adaptor{loc, title}); + return std::make_unique(loc, title); }); diff --git a/src/atlas/util/Bitflags.h b/src/atlas/util/Bitflags.h index dd0311663..b6e4f4f49 100644 --- a/src/atlas/util/Bitflags.h +++ b/src/atlas/util/Bitflags.h @@ -17,14 +17,6 @@ namespace util { //---------------------------------------------------------------------------------------------------------------------- -// Forward declaration of BitflagsView, defined further down in this file -namespace detail { -template -class BitflagsView {}; -} // namespace detail - -//---------------------------------------------------------------------------------------------------------------------- - /// @brief Convenience class to modify and interpret bitflags class Bitflags { public: @@ -53,12 +45,12 @@ class Bitflags { } /// @brief Create convenience accessor to modify flags - /// @note Use `auto` for return type! (will be more clear with C++14) - static detail::BitflagsView view(int& flags); + /// @note `auto` for return type! BitflagsView is a detail + static auto view(int& flags); /// @brief Create convenience accessor to modify flags - /// @note Use `auto` for return type! (will be more clear with C++14) - static detail::BitflagsView view(const int& flags); + /// @note `auto` for return type! BitflagsView is a detail + static auto view(const int& flags); }; //---------------------------------------------------------------------------------------------------------------------- @@ -67,6 +59,11 @@ namespace detail { //---------------------------------------------------------------------------------------------------------------------- +template +class BitflagsView {}; + +//---------------------------------------------------------------------------------------------------------------------- + // Template specialication for constant flags. There are no functions to edit the flags template <> class BitflagsView { @@ -103,10 +100,10 @@ class BitflagsView { //---------------------------------------------------------------------------------------------------------------------- -inline detail::BitflagsView Bitflags::view(const int& flags) { +inline auto Bitflags::view(const int& flags) { return detail::BitflagsView(flags); } -inline detail::BitflagsView Bitflags::view(int& flags) { +inline auto Bitflags::view(int& flags) { return detail::BitflagsView(flags); } From 81c4c2476478a9fa106848d61ab20f3d5c94a6f5 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 2 Mar 2023 15:55:53 +0000 Subject: [PATCH 03/19] Add missing prefix to test --- src/tests/functionspace/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/functionspace/CMakeLists.txt b/src/tests/functionspace/CMakeLists.txt index 318a7f285..47a1fff31 100644 --- a/src/tests/functionspace/CMakeLists.txt +++ b/src/tests/functionspace/CMakeLists.txt @@ -53,7 +53,7 @@ ecbuild_add_test( TARGET atlas_test_cubedsphere_functionspace ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -ecbuild_add_test( TARGET test_structuredcolumns_biperiodic +ecbuild_add_test( TARGET atlas_test_structuredcolumns_biperiodic SOURCES test_structuredcolumns_biperiodic.cc LIBS atlas MPI 5 From bcaeeeb85731ffa699c024216300fb120a819afb Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 2 Mar 2023 15:56:56 +0000 Subject: [PATCH 04/19] Workaround compiler bugs in aocc/4.0.0 --- atlas_io/tests/test_io_encoding.cc | 17 +++++++---------- .../partitioner/CheckerboardPartitioner.cc | 5 +++-- src/tests/io/test_io_encoding.cc | 15 +++++++-------- 3 files changed, 17 insertions(+), 20 deletions(-) diff --git a/atlas_io/tests/test_io_encoding.cc b/atlas_io/tests/test_io_encoding.cc index 9f327dc77..9cfccd52f 100644 --- a/atlas_io/tests/test_io_encoding.cc +++ b/atlas_io/tests/test_io_encoding.cc @@ -499,7 +499,6 @@ struct EncodedArray { std::vector in; }; - // ------------------------------------------------------------------------------------------------------- CASE("Decoding to std::vector") { @@ -570,18 +569,16 @@ CASE("Encode/Decode byte array") { auto validate = [&]() { EXPECT(out == encoded); - auto str = [](std::byte byte) { - std::bitset<8> bitset(reinterpret_cast(byte)); - return bitset.to_string(); + auto to_byte = []( const char* str) { + return std::byte(std::bitset<8>(str).to_ulong()); }; - EXPECT_EQ(str(out[0]), "00000001"); - EXPECT_EQ(str(out[1]), "00000011"); - EXPECT_EQ(str(out[2]), "00000111"); - EXPECT_EQ(str(out[3]), "00001111"); - EXPECT_EQ(str(out[4]), "00011111"); + EXPECT(out[0] == to_byte("00000001")); + EXPECT(out[1] == to_byte("00000011")); + EXPECT(out[2] == to_byte("00000111")); + EXPECT(out[3] == to_byte("00001111")); + EXPECT(out[4] == to_byte("00011111")); }; - SECTION("decode directly") { EXPECT_NO_THROW(decode(encoded.metadata, encoded.data, out)); validate(); diff --git a/src/atlas/grid/detail/partitioner/CheckerboardPartitioner.cc b/src/atlas/grid/detail/partitioner/CheckerboardPartitioner.cc index a06788be5..d452a99a2 100644 --- a/src/atlas/grid/detail/partitioner/CheckerboardPartitioner.cc +++ b/src/atlas/grid/detail/partitioner/CheckerboardPartitioner.cc @@ -63,8 +63,9 @@ CheckerboardPartitioner::Checkerboard CheckerboardPartitioner::checkerboard(cons } else { // default number of bands - double zz = std::sqrt((double)(nparts * cb.ny) / cb.nx); // aim at +/-square regions - cb.nbands = (idx_t)zz + 0.5; + double zz = std::sqrt(double(nparts * cb.ny) / double(cb.nx)); // aim at +/-square regions + cb.nbands = std::floor(zz + 0.5); + if (cb.nbands < 1) { cb.nbands = 1; // at least one band } diff --git a/src/tests/io/test_io_encoding.cc b/src/tests/io/test_io_encoding.cc index 2f9137b5e..7efb701b5 100644 --- a/src/tests/io/test_io_encoding.cc +++ b/src/tests/io/test_io_encoding.cc @@ -825,15 +825,14 @@ CASE("Encode/Decode byte array") { auto validate = [&]() { EXPECT(out == encoded); - auto str = [](std::byte byte) { - std::bitset<8> bitset(reinterpret_cast(byte)); - return bitset.to_string(); + auto to_byte = []( const char* str) { + return std::byte(std::bitset<8>(str).to_ulong()); }; - EXPECT_EQ(str(out[0]), "00000001"); - EXPECT_EQ(str(out[1]), "00000011"); - EXPECT_EQ(str(out[2]), "00000111"); - EXPECT_EQ(str(out[3]), "00001111"); - EXPECT_EQ(str(out[4]), "00011111"); + EXPECT(out[0] == to_byte("00000001")); + EXPECT(out[1] == to_byte("00000011")); + EXPECT(out[2] == to_byte("00000111")); + EXPECT(out[3] == to_byte("00001111")); + EXPECT(out[4] == to_byte("00011111")); }; From 60ccacca82d507d03cb81b24d7050383285b5beb Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 2 Mar 2023 16:27:14 +0000 Subject: [PATCH 05/19] Fix modernisation warning --- src/atlas/interpolation/method/cubedsphere/CellFinder.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/atlas/interpolation/method/cubedsphere/CellFinder.cc b/src/atlas/interpolation/method/cubedsphere/CellFinder.cc index c4c1e4d09..17f67f7d9 100644 --- a/src/atlas/interpolation/method/cubedsphere/CellFinder.cc +++ b/src/atlas/interpolation/method/cubedsphere/CellFinder.cc @@ -37,8 +37,8 @@ CellFinder::CellFinder(const Mesh& mesh, const util::Config& config): mesh_{mesh auto halo = config.getInt("halo", 0); for (idx_t i = 0; i < mesh_.cells().size(); ++i) { if (haloView(i) <= halo) { - points.push_back(PointLonLat(lonlatView(i, LON), lonlatView(i, LAT))); - payloads.push_back(i); + points.emplace_back(PointLonLat(lonlatView(i, LON), lonlatView(i, LAT))); + payloads.emplace_back(i); } } From 488a80d0735bf62c6fdc4c0c0ac34d371f86b943 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Du=C5=A1an=20Figala?= <34573508+figi44@users.noreply.github.com> Date: Tue, 7 Mar 2023 15:18:47 +0100 Subject: [PATCH 06/19] Enable CI on self-hosted runners (#117) For now this is without fckit and MPI --- .github/workflows/ci.yml | 75 +++++++++++++++++++++++++++ .github/workflows/label-public-pr.yml | 10 ++++ .github/workflows/reusable-ci-hpc.yml | 28 ++++++++++ .github/workflows/reusable-ci.yml | 28 ++++++++++ 4 files changed, 141 insertions(+) create mode 100644 .github/workflows/ci.yml create mode 100644 .github/workflows/label-public-pr.yml create mode 100644 .github/workflows/reusable-ci-hpc.yml create mode 100644 .github/workflows/reusable-ci.yml diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 000000000..28bd7a5c4 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,75 @@ +name: ci + +on: + # Trigger the workflow on push to master or develop, except tag creation + push: + branches: + - 'master' + - 'develop' + tags-ignore: + - '**' + + # Trigger the workflow on pull request + pull_request: ~ + + # Trigger the workflow manually + workflow_dispatch: ~ + + # Trigger after public PR approved for CI + pull_request_target: + types: [labeled] + +jobs: + # Run CI including downstream packages on self-hosted runners + downstream-ci: + name: downstream-ci + if: ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }} + uses: ecmwf-actions/downstream-ci/.github/workflows/downstream-ci.yml@main + with: + atlas: ecmwf/atlas@${{ github.event.pull_request.head.sha || github.sha }} + secrets: inherit + + # Run CI of private downstream packages on self-hosted runners + private-downstream-ci: + name: private-downstream-ci + needs: [downstream-ci] + if: (success() || failure()) && ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }} + runs-on: ubuntu-latest + permissions: + pull-requests: write + steps: + - name: Dispatch private downstream CI + uses: ecmwf-actions/dispatch-private-downstream-ci@v1 + with: + token: ${{ secrets.GH_REPO_READ_TOKEN }} + owner: ecmwf-actions + repository: private-downstream-ci + event_type: downstream-ci + payload: '{"atlas": "ecmwf/atlas@${{ github.event.pull_request.head.sha || github.sha }}"}' + + # Build downstream packages on HPC + downstream-ci-hpc: + name: downstream-ci-hpc + if: ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }} + uses: ecmwf-actions/downstream-ci/.github/workflows/downstream-ci-hpc.yml@main + with: + atlas: ecmwf/atlas@${{ github.event.pull_request.head.sha || github.sha }} + secrets: inherit + + # Run CI of private downstream packages on HPC + private-downstream-ci-hpc: + name: private-downstream-ci-hpc + needs: [downstream-ci-hpc] + if: (success() || failure()) && ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }} + runs-on: ubuntu-latest + permissions: + pull-requests: write + steps: + - name: Dispatch private downstream CI + uses: ecmwf-actions/dispatch-private-downstream-ci@v1 + with: + token: ${{ secrets.GH_REPO_READ_TOKEN }} + owner: ecmwf-actions + repository: private-downstream-ci + event_type: downstream-ci-hpc + payload: '{"atlas": "ecmwf/atlas@${{ github.event.pull_request.head.sha || github.sha }}"}' diff --git a/.github/workflows/label-public-pr.yml b/.github/workflows/label-public-pr.yml new file mode 100644 index 000000000..59b2bfa2b --- /dev/null +++ b/.github/workflows/label-public-pr.yml @@ -0,0 +1,10 @@ +# Manage labels of pull requests that originate from forks +name: label-public-pr + +on: + pull_request_target: + types: [opened, synchronize] + +jobs: + label: + uses: ecmwf-actions/reusable-workflows/.github/workflows/label-pr.yml@v2 diff --git a/.github/workflows/reusable-ci-hpc.yml b/.github/workflows/reusable-ci-hpc.yml new file mode 100644 index 000000000..5f4c5d15f --- /dev/null +++ b/.github/workflows/reusable-ci-hpc.yml @@ -0,0 +1,28 @@ +name: reusable-ci-hpc + +on: + workflow_call: + inputs: + atlas: + required: false + type: string + eckit: + required: false + type: string +jobs: + ci-hpc: + name: ci-hpc + uses: ecmwf-actions/reusable-workflows/.github/workflows/ci-hpc.yml@v2 + with: + name-prefix: atlas- + build-inputs: | + --package: ${{ inputs.atlas || 'ecmwf/atlas@develop' }} + --modules: | + ecbuild + ninja + --modules-package: | + openmpi + --dependencies: | + ${{ inputs.eckit || 'ecmwf/eckit@develop' }} + --parallel: 64 + secrets: inherit diff --git a/.github/workflows/reusable-ci.yml b/.github/workflows/reusable-ci.yml new file mode 100644 index 000000000..e6a5aac35 --- /dev/null +++ b/.github/workflows/reusable-ci.yml @@ -0,0 +1,28 @@ +name: reusable-ci + +on: + workflow_call: + inputs: + atlas: + required: false + type: string + eckit: + required: false + type: string + +jobs: + ci: + name: atlas-ci + uses: ecmwf-actions/reusable-workflows/.github/workflows/ci.yml@v2 + with: + repository: ${{ inputs.atlas || 'ecmwf/atlas@develop' }} + name_prefix: atlas- + build_package_inputs: | + repository: ${{ inputs.atlas || 'ecmwf/atlas@develop' }} + self_coverage: true + dependencies: | + ecmwf/ecbuild + ${{ inputs.eckit || 'ecmwf/eckit' }} + dependency_branch: develop + parallelism_factor: 8 + secrets: inherit From ac5ba8ca909db226a564ba49b543d437ae08c123 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Du=C5=A1an=20Figala?= <34573508+figi44@users.noreply.github.com> Date: Wed, 8 Mar 2023 15:16:19 +0100 Subject: [PATCH 07/19] Fix `More processors requested than permitted` error on HPC (#118) * Add fftw and eigen modules * Remove openmpi --- .github/workflows/reusable-ci-hpc.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/reusable-ci-hpc.yml b/.github/workflows/reusable-ci-hpc.yml index 5f4c5d15f..744456177 100644 --- a/.github/workflows/reusable-ci-hpc.yml +++ b/.github/workflows/reusable-ci-hpc.yml @@ -21,7 +21,8 @@ jobs: ecbuild ninja --modules-package: | - openmpi + fftw + eigen --dependencies: | ${{ inputs.eckit || 'ecmwf/eckit@develop' }} --parallel: 64 From 6d6b6e1adaba9c54317c03b92aedf2b2244cc1c4 Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Mon, 13 Mar 2023 15:00:59 +0000 Subject: [PATCH 08/19] atlas --info shows if built with PROJ --- src/CMakeLists.txt | 6 ++++++ src/atlas/library/Library.cc | 2 ++ src/atlas/library/defines.h.in | 1 + 3 files changed, 9 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 74c80f31f..4199ca74a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,6 +24,12 @@ else() set( atlas_HAVE_TESSELATION 0 ) endif() +if( atlas_HAVE_PROJ ) + set( atlas_HAVE_PROJ 1 ) +else() + set( atlas_HAVE_PROJ 0 ) +endif() + if( atlas_HAVE_FORTRAN ) set( atlas_HAVE_FORTRAN 1 ) else() diff --git a/src/atlas/library/Library.cc b/src/atlas/library/Library.cc index 0ce439bb6..66621621e 100644 --- a/src/atlas/library/Library.cc +++ b/src/atlas/library/Library.cc @@ -454,6 +454,7 @@ void Library::Information::print(std::ostream& out) const { bool feature_FFTW(ATLAS_HAVE_FFTW); bool feature_Eigen(ATLAS_HAVE_EIGEN); bool feature_Tesselation(ATLAS_HAVE_TESSELATION); + bool feature_PROJ(ATLAS_HAVE_PROJ); bool feature_BoundsChecking(ATLAS_ARRAYVIEW_BOUNDS_CHECKING); bool feature_Init_sNaN(ATLAS_INIT_SNAN); bool feature_MPI(false); @@ -478,6 +479,7 @@ void Library::Information::print(std::ostream& out) const { << " Eigen : " << str(feature_Eigen) << '\n' << " MKL : " << str(feature_MKL()) << '\n' << " Tesselation : " << str(feature_Tesselation) << '\n' + << " PROJ : " << str(feature_PROJ) << '\n' << " ArrayDataStore : " << array_data_store << '\n' << " idx_t : " << ATLAS_BITS_LOCAL << " bit integer" << '\n' << " gidx_t : " << ATLAS_BITS_GLOBAL << " bit integer" << '\n'; diff --git a/src/atlas/library/defines.h.in b/src/atlas/library/defines.h.in index f59630451..5056fafc8 100644 --- a/src/atlas/library/defines.h.in +++ b/src/atlas/library/defines.h.in @@ -23,6 +23,7 @@ #define ATLAS_HAVE_EIGEN @atlas_HAVE_EIGEN@ #define ATLAS_HAVE_FFTW @atlas_HAVE_FFTW@ #define ATLAS_HAVE_MPI @atlas_HAVE_MPI@ +#define ATLAS_HAVE_PROJ @atlas_HAVE_PROJ@ #define ATLAS_BITS_GLOBAL @ATLAS_BITS_GLOBAL@ #define ATLAS_ARRAYVIEW_BOUNDS_CHECKING @atlas_HAVE_BOUNDSCHECKING@ #define ATLAS_INDEXVIEW_BOUNDS_CHECKING @atlas_HAVE_BOUNDSCHECKING@ From b9efb25ca29d949fbcfb8f29328d9a5a5ebb6d25 Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Mon, 13 Mar 2023 15:01:18 +0000 Subject: [PATCH 09/19] atlas --info removal of unnecessary header --- src/atlas/library/Library.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/atlas/library/Library.cc b/src/atlas/library/Library.cc index 66621621e..86c7f1075 100644 --- a/src/atlas/library/Library.cc +++ b/src/atlas/library/Library.cc @@ -58,7 +58,6 @@ static bool feature_MKL() { #include "atlas/util/Config.h" #if ATLAS_HAVE_TRANS -#include "atlas/library/config.h" #if ATLAS_HAVE_ECTRANS #include "ectrans/transi.h" #else From 0bbb4e57c5127e38f739240c55ac9516fee73abe Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 24 Mar 2023 10:44:24 +0100 Subject: [PATCH 10/19] Avoid failing atlas_test_interpolation_structured2D_to_unstructured with more than 2 MPI tasks --- ...erpolation_structured2D_to_unstructured.cc | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc b/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc index 65ac5e42b..e3d99cf30 100644 --- a/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc +++ b/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc @@ -116,6 +116,9 @@ FieldSet create_target_fields(FunctionSpace& fs, idx_t nb_fields, idx_t nb_level CASE("test_match") { + if(mpi::size() > 2) { + return; + } idx_t nb_fields = 2; idx_t nb_levels = 3; @@ -132,6 +135,10 @@ CASE("test_match") { interpolation.execute(fields_source, fields_target); } CASE("test_nomatch") { + if(mpi::size() > 2) { + return; + } + idx_t nb_fields = 2; idx_t nb_levels = 3; @@ -151,6 +158,29 @@ CASE("test_nomatch") { } } +CASE("test_nomatch 2") { + if(mpi::size() > 2) { + return; + } + + idx_t nb_fields = 2; + idx_t nb_levels = 3; + + Grid input_grid(input_gridname("O32")); + StructuredColumns input_fs(input_grid, option::halo(1) | option::levels(nb_levels)); + + FunctionSpace output_fs = output_functionspace_nomatch(); + + if (false) // expected to throw + { + Interpolation interpolation(option::type("structured-linear2D"), input_fs, output_fs); + + FieldSet fields_source = create_source_fields(input_fs, nb_fields, nb_levels); + FieldSet fields_target = create_target_fields(output_fs, nb_fields, nb_levels); + + interpolation.execute(fields_source, fields_target); + } +} } // namespace test } // namespace atlas From cb9b213ecd01dedfb2529352d809754e876f09d8 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 27 Mar 2023 15:45:59 +0200 Subject: [PATCH 11/19] Support formatting for Config::json() --- src/atlas/util/Config.cc | 4 ++-- src/atlas/util/Config.h | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/atlas/util/Config.cc b/src/atlas/util/Config.cc index 64955b421..8ccf234c1 100644 --- a/src/atlas/util/Config.cc +++ b/src/atlas/util/Config.cc @@ -109,9 +109,9 @@ std::vector Config::keys() const { return result; } -std::string Config::json() const { +std::string Config::json(eckit::JSON::Formatting formatting) const { std::stringstream json; - eckit::JSON js(json); + eckit::JSON js(json,formatting); js << *this; return json.str(); } diff --git a/src/atlas/util/Config.h b/src/atlas/util/Config.h index 02170b307..97a736872 100644 --- a/src/atlas/util/Config.h +++ b/src/atlas/util/Config.h @@ -13,6 +13,7 @@ #include #include "eckit/config/LocalConfiguration.h" +#include "eckit/log/JSON.h" namespace eckit { class PathName; @@ -77,7 +78,7 @@ class Config : public eckit::LocalConfiguration { std::vector keys() const; - std::string json() const; + std::string json(eckit::JSON::Formatting = eckit::JSON::Formatting::indent()) const; }; // ------------------------------------------------------------------ From f07b1c3f6ddb5db3e9e31b538b18461136d52e53 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 23 Mar 2023 14:02:19 +0100 Subject: [PATCH 12/19] Add support for StructuredPartitionPolygon with halo > 0 --- .../functionspace/detail/StructuredColumns.cc | 19 ++- .../functionspace/detail/StructuredColumns.h | 4 +- src/atlas/grid/StructuredPartitionPolygon.cc | 132 +++++++++++------- 3 files changed, 93 insertions(+), 62 deletions(-) diff --git a/src/atlas/functionspace/detail/StructuredColumns.cc b/src/atlas/functionspace/detail/StructuredColumns.cc index 140e47c0f..02bd6cbc9 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.cc +++ b/src/atlas/functionspace/detail/StructuredColumns.cc @@ -485,18 +485,23 @@ void StructuredColumns::Map2to1::resize(std::array i_range, std::array } const util::PartitionPolygon& StructuredColumns::polygon(idx_t halo) const { - if (not polygon_) { - polygon_.reset(new grid::StructuredPartitionPolygon(*this, halo)); + if (halo >= static_cast(polygons_.size())) { + polygons_.resize(halo + 1); } - return *polygon_; + if (not polygons_[halo]) { + if (halo > this->halo()) { + throw_Exception("StructuredColumsn does not contain a halo of size " + std::to_string(halo) + ".", Here()); + } + polygons_[halo].reset(new grid::StructuredPartitionPolygon(*this, halo)); + } + return *polygons_[halo]; } const util::PartitionPolygons& StructuredColumns::polygons() const { - if (polygons_.size()) { - return polygons_; + if (all_polygons_.size() == 0) { + polygon().allGather(all_polygons_); } - polygon().allGather(polygons_); - return polygons_; + return all_polygons_; } // ---------------------------------------------------------------------------- diff --git a/src/atlas/functionspace/detail/StructuredColumns.h b/src/atlas/functionspace/detail/StructuredColumns.h index f3d3d6850..392ed577c 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.h +++ b/src/atlas/functionspace/detail/StructuredColumns.h @@ -218,8 +218,8 @@ class StructuredColumns : public FunctionSpaceImpl { mutable util::ObjectHandle gather_scatter_; mutable util::ObjectHandle checksum_; mutable util::ObjectHandle halo_exchange_; - mutable std::unique_ptr polygon_; - mutable util::PartitionPolygons polygons_; + mutable std::vector> polygons_; + mutable util::PartitionPolygons all_polygons_; Field field_xy_; Field field_partition_; diff --git a/src/atlas/grid/StructuredPartitionPolygon.cc b/src/atlas/grid/StructuredPartitionPolygon.cc index e5e10b816..e5b3a6773 100644 --- a/src/atlas/grid/StructuredPartitionPolygon.cc +++ b/src/atlas/grid/StructuredPartitionPolygon.cc @@ -35,10 +35,32 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect const auto grid = fs.grid(); const auto dom = RectangularDomain(grid.domain()); - if (_halo > 0) { - throw_Exception("halo must be smaller than that of StructuredColumns", Here()); + if (_halo > 0 && _halo != fs.halo()) { + throw_Exception("halo must zero or matching the one of the StructuredColumns", Here()); } + bool north_pole_included = 90. - grid.y(0) == 0.; + bool south_pole_included = 90. + grid.y(grid.ny() - 1) == 0; + + auto compute_j = [&](const idx_t j) { + idx_t jj; + if (j < 0) { + return -j - 1 + north_pole_included; + } + else if (j >= grid.ny()) { + return 2 * grid.ny() - j - 1 - south_pole_included; + } + return j; + }; + + auto compute_y = [&](const idx_t j) { + return fs.compute_xy(0,j).y(); + }; + + auto compute_x = [&](const idx_t i, const idx_t j) { + return fs.compute_xy(i,j).x(); + }; + auto equal = [](const double& a, const double& b) { return std::abs(a - b) < 1.e-12; }; auto last_edge_horizontal = [&]() { @@ -58,11 +80,10 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect PointXY p; idx_t c{0}; - idx_t i, j; auto add_point = [&](const Point2& point) { points.emplace_back(point); - // Log::info() << "add point (" << points.size() - 1 << ") " << point << std::endl; + // Log::info() << "add point (" << points.size() - 1 << ") " << point << std::endl; }; auto add_vertical_edge = [&](const Point2& point) { @@ -91,38 +112,38 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect // Top // Top left point { - j = fs.j_begin(); - i = fs.i_begin(j); - if (j == 0) { + idx_t j = _halo ? fs.j_begin_halo() : fs.j_begin(); + idx_t i = _halo ? fs.i_begin_halo(j) : fs.i_begin(j); + if (j == 0 && _halo == 0) { p[YY] = dom.ymax(); } else { - p[YY] = 0.5 * (grid.y(j - 1) + grid.y(j)); + p[YY] = 0.5 * (compute_y(j - 1) + compute_y(j)); } - if (i == 0) { + if (i == 0 && _halo == 0) { p[XX] = dom.xmin(); } else { - p[XX] = 0.5 * (grid.x(i - 1, j) + grid.x(i, j)); + p[XX] = 0.5 * (compute_x(i - 1, j) + compute_x(i, j)); } add_point(p); } // Top right point { - j = fs.j_begin(); - i = fs.i_end(j) - 1; - if (j == 0) { + idx_t j = _halo ? fs.j_begin_halo() : fs.j_begin(); + idx_t i = _halo ? fs.i_end_halo(j) - 1 : fs.i_end(j) - 1; + if (j == 0 && _halo == 0) { p[YY] = dom.ymax(); } else { - p[YY] = 0.5 * (grid.y(j - 1) + grid.y(j)); + p[YY] = 0.5 * (compute_y(j - 1) + compute_y(j)); } - if (i == grid.nx(j) - 1) { + if (i == grid.nx(compute_j(j)) - 1 && _halo == 0) { p[XX] = dom.xmax(); } else { - p[XX] = 0.5 * (grid.x(i, j) + grid.x(i + 1, j)); + p[XX] = 0.5 * (compute_x(i, j) + compute_x(i + 1, j)); } add_horizontal_edge(p); @@ -131,14 +152,16 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect } // Right side { - for (j = fs.j_begin(); j < fs.j_end() - 1; ++j) { - p[YY] = grid.y(j); - i = fs.i_end(j) - 1; - if (i == grid.nx(j) - 1) { + idx_t j_begin = _halo ? fs.j_begin_halo() : fs.j_begin(); + idx_t j_end = _halo ? fs.j_end_halo() : fs.j_end(); + for (idx_t j = j_begin; j < j_end - 1 ; ++j) { + p[YY] = compute_y(j); + idx_t i = _halo ? fs.i_end_halo(j) - 1 : fs.i_end(j) - 1; + if (i == grid.nx(compute_j(j)) - 1 && _halo == 0) { p[XX] = dom.xmax(); } else { - p[XX] = 0.5 * (grid.x(i, j) + grid.x(i + 1, j)); + p[XX] = 0.5 * (compute_x(i, j) + compute_x(i + 1, j)); } if (p == points.back()) { continue; @@ -168,19 +191,20 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect // Bottom // Bottom right point(s) { - j = fs.j_end() - 1; - i = fs.i_end(j) - 1; - if (j == grid.ny() - 1) { + idx_t j = _halo ? fs.j_end_halo() - 1 : fs.j_end() - 1; + idx_t i = _halo ? fs.i_end_halo(j) - 1 : fs.i_end(j) - 1; + + if (j == grid.ny() - 1 && _halo == 0) { p[YY] = dom.ymin(); } else { - p[YY] = 0.5 * (grid.y(j) + grid.y(j + 1)); + p[YY] = 0.5 * (compute_y(j) + compute_y(j + 1)); } - if (i == grid.nx(j) - 1) { + if (i == grid.nx(compute_j(j)) - 1 && _halo == 0) { p[XX] = dom.xmax(); } else { - p[XX] = 0.5 * (grid.x(i, j) + grid.x(i + 1, j)); + p[XX] = 0.5 * (compute_x(i, j) + compute_x(i + 1, j)); } ymin = p[YY]; @@ -191,7 +215,7 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect PointXY ptmp; ptmp = p; p[XX] = points.back()[XX]; - p[YY] = 0.5 * (points.back()[YY] + grid.y(j)); + p[YY] = 0.5 * (points.back()[YY] + compute_y(j)); add_vertical_edge(p); pmin = p; @@ -202,11 +226,12 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect ptmp = p; p[YY] = points.back()[YY]; + add_horizontal_edge(p); p = ptmp; } - if (xmax - grid.xspace().dx()[j] < grid.x(i, j)) { - xmax = std::min(xmax, 0.5 * (grid.x(i + 1, j) + grid.x(i, j))); + if (xmax - grid.xspace().dx()[compute_j(j)] < compute_x(i, j)) { + xmax = std::min(xmax, 0.5 * (compute_x(i + 1, j) + compute_x(i, j))); } else { ymin = pmin[YY]; @@ -216,46 +241,48 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect // Bottom left point { - j = fs.j_end() - 1; - i = fs.i_begin(j); - if (j == grid.ny() - 1) { + idx_t j = _halo ? fs.j_end_halo() - 1 : fs.j_end() - 1; + idx_t i = _halo ? fs.i_begin_halo(j) : fs.i_begin(j); + if (j == grid.ny() - 1 && _halo == 0) { p[YY] = dom.ymin(); } else { - p[YY] = 0.5 * (grid.y(j) + grid.y(j + 1)); + p[YY] = 0.5 * (compute_y(j) + compute_y(j + 1)); } - if (i == 0) { + if (i == 0 && _halo == 0) { p[XX] = dom.xmin(); } else { - p[XX] = 0.5 * (grid.x(i - 1, j) + grid.x(i, j)); + p[XX] = 0.5 * (compute_x(i - 1, j) + compute_x(i, j)); } add_horizontal_edge(p); xmin = p[XX]; } // Left side { - for (j = fs.j_end() - 1; j >= fs.j_begin(); --j) { - p[YY] = grid.y(j); - i = fs.i_begin(j); - if (i == 0) { + idx_t j_end = _halo ? fs.j_end_halo() : fs.j_end(); + idx_t j_begin = _halo ? fs.j_begin_halo() : fs.j_begin(); + for (idx_t j = j_end - 1; j >= j_begin; --j) { + p[YY] = compute_y(j); + idx_t i = _halo ? fs.i_begin_halo(j) : fs.i_begin(j); + if (i == 0 && _halo == 0) { p[XX] = dom.xmin(); } else { - p[XX] = 0.5 * (grid.x(i - 1, j) + grid.x(i, j)); + p[XX] = 0.5 * (compute_x(i - 1, j) + compute_x(i, j)); } - if (j > fs.j_begin()) { + if (j > j_begin) { xmin = std::max(xmin, p[XX]); } - if (j == fs.j_begin() + 1 && xmin < points[0][XX]) { - idx_t jtop = fs.j_begin(); - idx_t itop = fs.i_begin(jtop); - if (xmin + grid.xspace().dx()[jtop] > grid.x(itop, jtop)) { - xmin = std::max(xmin, 0.5 * (grid.x(itop - 1, jtop) + grid.x(itop, jtop))); + if (j == j_begin && xmin < points[0][XX]) { + idx_t jtop = j_begin; + idx_t itop = _halo ? fs.i_begin_halo(jtop) : fs.i_begin(jtop); + if (xmin + grid.xspace().dx()[compute_j(jtop)] > compute_x(itop, jtop)) { + xmin = std::max(xmin, 0.5 * (compute_x(itop - 1, jtop) + compute_x(itop, jtop))); } else { - ymax = 0.5 * (p[YY] + grid.y(jtop)); + ymax = 0.5 * (p[YY] + compute_y(jtop)); } } @@ -269,7 +296,7 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect // vertical edge p[XX] = points.back()[XX]; - p[YY] = 0.5 * (points.back()[YY] + grid.y(j)); + p[YY] = 0.5 * (points.back()[YY] + compute_y(j)); add_vertical_edge(p); p = ptmp; @@ -358,8 +385,7 @@ void StructuredPartitionPolygon::outputPythonScript(const eckit::PathName& filep "\n" "from itertools import cycle" "\n" "import matplotlib.cm as cm" "\n" "import numpy as np" - "\n" "cycol = cycle([cm.Paired(i) for i in " - "np.linspace(0,1,12,endpoint=True)]).next" + "\n" "colours = cycle([cm.Paired(i) for i in np.linspace(0,1,12,endpoint=True)])" "\n" "\n" "fig = plt.figure()" "\n" "ax = fig.add_subplot(111,aspect='equal')" @@ -391,8 +417,8 @@ void StructuredPartitionPolygon::outputPythonScript(const eckit::PathName& filep f << "]"; } f << "\n" - "\n" "c = cycol()" - "\n" "ax.add_patch(patches.PathPatch(Path(verts_" << r << ", codes_" << r << "), facecolor=c, color=c, alpha=0.3, lw=1))"; + "\n" "c = next(colours)" + "\n" "ax.add_patch(patches.PathPatch(Path(verts_" << r << ", codes_" << r << "), edgecolor=c, facecolor=c, alpha=0.3, lw=1))"; if ( plot_nodes ) { f << "\n" "if plot_nodes:" "\n" " ax.scatter(x_" << r << ", y_" << r << ", color=c, marker='o')"; From 929d65e48f29e3349795ac0a4896a511327af606 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 23 Mar 2023 14:03:11 +0100 Subject: [PATCH 13/19] Add PartitionPolygon::json() function --- src/atlas/util/Polygon.cc | 30 ++++++++++++++++++++++++++++++ src/atlas/util/Polygon.h | 8 +++++++- 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/atlas/util/Polygon.cc b/src/atlas/util/Polygon.cc index 40a88cb63..dc5b38872 100644 --- a/src/atlas/util/Polygon.cc +++ b/src/atlas/util/Polygon.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include "eckit/types/FloatCompare.h" @@ -280,6 +281,35 @@ void ExplicitPartitionPolygon::allGather(PartitionPolygons&) const { ATLAS_NOTIMPLEMENTED; } +std::string PartitionPolygons::json(const eckit::Configuration& config) const { + std::ostringstream out; + out << "[\n"; + for (size_t i=0; i {}; +class PartitionPolygons : public VectorOfAbstract { +public: + std::string json(const eckit::Configuration& = util::NoConfig()) const; +}; //------------------------------------------------------------------------------------------------------ From b2b482c9a828ea438842316433d1b474f9a66872 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 24 Mar 2023 10:45:45 +0100 Subject: [PATCH 14/19] Add test for StructuredInterpolation2D to show normalisation issue --- src/tests/interpolation/CMakeLists.txt | 10 ++ ...st_interpolation_structured2D_to_points.cc | 150 ++++++++++++++++++ 2 files changed, 160 insertions(+) create mode 100644 src/tests/interpolation/test_interpolation_structured2D_to_points.cc diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index f2cbda9a2..24eb949cc 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -102,6 +102,16 @@ ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_unstructured ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) + +ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_points + SOURCES test_interpolation_structured2D_to_points.cc + LIBS atlas + MPI 4 + CONDITION eckit_HAVE_MPI + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + + ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere SOURCES test_interpolation_cubedsphere.cc LIBS atlas diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_points.cc b/src/tests/interpolation/test_interpolation_structured2D_to_points.cc new file mode 100644 index 000000000..956520ee0 --- /dev/null +++ b/src/tests/interpolation/test_interpolation_structured2D_to_points.cc @@ -0,0 +1,150 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/PointCloud.h" +#include "atlas/functionspace/StructuredColumns.h" +#include "atlas/grid/Grid.h" +#include "atlas/grid/Iterator.h" +#include "atlas/interpolation.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" +#include "atlas/util/PolygonXY.h" + +#include "tests/AtlasTestEnvironment.h" + +using atlas::functionspace::PointCloud; +using atlas::functionspace::StructuredColumns; +using atlas::util::Config; + +namespace atlas { +namespace test { + +//----------------------------------------------------------------------------- + +std::string input_gridname(const std::string& default_grid) { + return eckit::Resource("--input-grid", default_grid); +} + +FunctionSpace output_functionspace( const FunctionSpace& input_functionspace, bool expect_fail ) { + + std::vector all_points { + {360., 90.}, + {360., 0.}, + {360., -90.}, + {0.,0.}, + }; + + // Only keep points that match the input partitioning. + // Note that with following implementation it could be that some points + // are present in two partitions, but it is not a problem for this test purpose. + std::vector points; + auto polygon = util::PolygonXY{input_functionspace.polygon()}; + for (const auto& p : all_points ) { + if (polygon.contains(p)) { + points.emplace_back(p); + } + } + if( expect_fail && mpi::rank() == mpi::size() - 1 ) { + points.emplace_back(720,0.); + } + return PointCloud(points); +} + + +FieldSet create_source_fields(StructuredColumns& fs, idx_t nb_fields, idx_t nb_levels) { + using Value = double; + FieldSet fields_source; + auto lonlat = array::make_view(fs.xy()); + for (idx_t f = 0; f < nb_fields; ++f) { + auto field_source = fields_source.add(fs.createField()); + auto source = array::make_view(field_source); + for (idx_t n = 0; n < fs.size(); ++n) { + for (idx_t k = 0; k < nb_levels; ++k) { + source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + } + }; + } + return fields_source; +} +FieldSet create_target_fields(FunctionSpace& fs, idx_t nb_fields, idx_t nb_levels) { + using Value = double; + FieldSet fields_target; + for (idx_t f = 0; f < nb_fields; ++f) { + fields_target.add(fs.createField(option::levels(nb_levels))); + } + return fields_target; +} + +void do_test( std::string type, int input_halo, bool matrix_free, bool expect_failure ) { + idx_t nb_fields = 2; + idx_t nb_levels = 3; + + Grid input_grid(input_gridname("O32")); + StructuredColumns input_fs(input_grid, option::levels(nb_levels) | + option::halo(input_halo)); + + FunctionSpace output_fs = output_functionspace(input_fs, expect_failure); + + Interpolation interpolation(option::type(type) | + util::Config("matrix_free",matrix_free) | + util::Config("verbose",eckit::Resource("--verbose",false)), + input_fs, output_fs); + + FieldSet fields_source = create_source_fields(input_fs, nb_fields, nb_levels); + FieldSet fields_target = create_target_fields(output_fs, nb_fields, nb_levels); + + interpolation.execute(fields_source, fields_target); +} + +CASE("test structured-bilinear, halo 2, with matrix") { + EXPECT_NO_THROW( do_test("structured-bilinear",2,false,false) ); +} + +CASE("test structured-bilinear, halo 2, without matrix") { + EXPECT_NO_THROW( do_test("structured-bilinear",2,true,false) ); +} + +CASE("test structured-bilinear, halo 2, with matrix, expected failure") { + EXPECT_THROWS_AS( do_test("structured-bilinear",2,false,true), eckit::Exception ); +} + +CASE("test structured-bilinear, halo 2, without matrix, expected failure") { + EXPECT_THROWS_AS( do_test("structured-bilinear",2,false,true), eckit::Exception ); +} + +CASE("test structured-bilinear, halo 1, with matrix, expected failure") { + EXPECT_THROWS_AS( do_test("structured-bilinear",1,false,false), eckit::Exception ); +} + +CASE("test structured-bilinear, halo 1, without matrix, expected failure") { + EXPECT_THROWS_AS( do_test("structured-bilinear",1,true,false), eckit::Exception ); +} + +CASE("test structured-bicubic, halo 3, with matrix") { + EXPECT_NO_THROW( do_test("structured-bicubic",3,false,false) ); +} + +CASE("test structured-bicubic, halo 2, with matrix") { + EXPECT_THROWS_AS( do_test("structured-bicubic",2,false,false), eckit::Exception ); +} + + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} From ea6400aa8358ea7814ff3234f84c51a03e29f61e Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 23 Mar 2023 14:13:20 +0100 Subject: [PATCH 15/19] Fix StructuredInterpolation2D with retry for failed stencils --- .../structured/StructuredInterpolation2D.h | 16 +- .../structured/StructuredInterpolation2D.tcc | 419 ++++++++++++++---- .../structured/StructuredInterpolation3D.tcc | 6 +- .../method/structured/kernels/Cubic3DKernel.h | 4 +- .../kernels/CubicHorizontalKernel.h | 41 +- .../structured/kernels/Linear3DKernel.h | 4 +- .../kernels/LinearHorizontalKernel.h | 43 +- .../structured/kernels/QuasiCubic3DKernel.h | 4 +- .../kernels/QuasiCubicHorizontalKernel.h | 41 +- .../test_interpolation_structured2D.cc | 1 + 10 files changed, 451 insertions(+), 128 deletions(-) diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h index 9cd9ffb0c..77598fefc 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h @@ -35,17 +35,18 @@ class StructuredInterpolation2D : public Method { public: StructuredInterpolation2D(const Config& config); - virtual ~StructuredInterpolation2D() override {} + ~StructuredInterpolation2D() override {} - virtual void print(std::ostream&) const override; + void print(std::ostream&) const override; + const Kernel& kernel() const { return *kernel_; } -protected: - void setup(const FunctionSpace& source); + const FunctionSpace& source() const override { return source_; } - virtual const FunctionSpace& source() const override { return source_; } + const FunctionSpace& target() const override { return target_; } - virtual const FunctionSpace& target() const override { return target_; } +protected: + void setup(const FunctionSpace& source); private: virtual void do_setup(const Grid& source, const Grid& target, const Cache&) override; @@ -75,6 +76,9 @@ class StructuredInterpolation2D : public Method { FunctionSpace target_; bool matrix_free_; + bool verbose_; + double convert_units_; + idx_t out_npts_; std::unique_ptr kernel_; }; diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index 4375ba4e0..8f9ae1392 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -12,6 +12,16 @@ #include "StructuredInterpolation2D.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include "eckit/filesystem/LocalPathName.h" #include "atlas/array/ArrayView.h" #include "atlas/field/Field.h" @@ -35,6 +45,178 @@ namespace atlas { namespace interpolation { namespace method { +namespace structured2d { + +namespace { + + inline std::string search_replace(const std::string& in, const std::string& search, const std::string& replace) { + std::string out = in; + int pos = out.find(search); + while (pos != std::string::npos) { + out.erase(pos, search.length()); + out.insert(pos, replace); + pos = out.find(search, pos + replace.length()); + } + return out; + } + + inline std::string left_padded_str(int x, int max = 0) { + auto digits = [](long x) -> long { return std::floor(std::log10(std::max(1l, x))) + 1l; }; + if (max) { + std::ostringstream ss; + ss << std::setw(digits(max)) << std::setfill('0') << x; + return ss.str(); + } + return std::to_string(x); + } + + inline std::string filename( const std::string& path ){ + return search_replace(path, "%p", left_padded_str(mpi::rank(),mpi::size())); + } + + template + inline std::string to_str (const std::vector& points, LonLat lonlat) { + std::ostringstream out; + out << "[\n"; + for( size_t j=0; j::max(); + double src_east = std::numeric_limits::min(); + for( auto& p : src_polygon.lonlat() ) { + src_west = std::min(src_west, p[0]); + src_east = std::max(src_east, p[0]); + } + return std::make_tuple(src_west,src_east); + }; + + + inline void output_source_partition_polygon(const std::string& path, functionspace::StructuredColumns src_fs, idx_t halo) { + auto& polygon = src_fs.polygon(halo); + std::ofstream f(filename(path)); + f << "[\n " << polygon.json() << "\n]"; + } + + inline void output_all_source_partition_polygons(const std::string& path, functionspace::StructuredColumns src_fs, idx_t halo) { + auto& polygon = src_fs.polygon(halo); + util::PartitionPolygons polygons; + polygon.allGather(polygons); + if( mpi::rank() == 0 ) { + std::ofstream f(path); + f << polygons.json(); + } + } + + inline void output_source_grid_points(const std::string& path, functionspace::StructuredColumns src_fs, idx_t halo) { + std::ofstream f(filename(path)); + f << "[\n"; + if( halo == 0 ) { + idx_t size = src_fs.sizeOwned(); + idx_t n = 0; + for( idx_t j=src_fs.j_begin(); j + inline void output_target_points(const std::string& path, const std::vector& points, LonLat lonlat) { + std::ofstream f(filename(path)); + f << to_str(points,lonlat); + } + + template + inline void handle_failed_points(const Interpolation& interpolation, const std::vector& failed_points, LonLat lonlat ) { + const auto src_fs = functionspace::StructuredColumns( interpolation.source() ); + size_t num_failed_points{0}; + mpi::comm().allReduce(failed_points.size(), num_failed_points, eckit::mpi::sum()); + if( num_failed_points > 0 ) { + std::string halo_str = std::to_string(src_fs.halo()); + output_source_partition_polygon("atlas_source_partition_polygons_halo_0_p%p.json",src_fs,0); + output_source_partition_polygon("atlas_source_partition_polygons_halo_" + halo_str + "_p%p.json",src_fs,src_fs.halo()); + output_target_points("atlas_target_failed_points_p%p.json", failed_points, lonlat); + auto [src_west, src_east] = compute_src_west_east(src_fs); + idx_t my_rank = mpi::rank(); + for( idx_t p=0; p"); + f << src_fs_config.json(); + } + std::string p_range = "{"+left_padded_str(0,mpi::size()-1)+".."+left_padded_str(mpi::size()-1,mpi::size()-1)+"}"; + Log::info() << "Dumped files to " << eckit::LocalPathName::cwd() << ": \n" + << " atlas_target_failed_points_p" << p_range << ".json\n" + << " atlas_source_partition_polygons_halo_0_p" << p_range << ".json\n" + << " atlas_source_partition_polygons_halo_0.json\n" + << " atlas_source_partition_polygons_halo_" << src_fs.halo() << "_p" << p_range << ".json\n" + << " atlas_source_partition_polygons_halo_" << src_fs.halo() << ".json\n" + << " atlas_interpolation_info.json\n" + << std::endl; + + std::ostringstream err; + err << "StructuredInterpolation2D<" << interpolation.kernel().className() << "> failed for " + << num_failed_points << " points with source halo="< double StructuredInterpolation2D::convert_units_multiplier( const Field& field ) { std::string units = field.metadata().getString( "units", "degrees" ); @@ -50,8 +232,10 @@ double StructuredInterpolation2D::convert_units_multiplier( const Field& template StructuredInterpolation2D::StructuredInterpolation2D( const Method::Config& config ) : Method( config ), - matrix_free_{false} { + matrix_free_{false} , + verbose_{false} { config.get( "matrix_free", matrix_free_ ); + config.get( "verbose", verbose_ ); } @@ -140,61 +324,107 @@ void StructuredInterpolation2D::print( std::ostream& ) const { template void StructuredInterpolation2D::setup( const FunctionSpace& source ) { + using namespace structured2d; + kernel_.reset( new Kernel( source ) ); if ( functionspace::StructuredColumns( source ).halo() < 1 ) { throw_Exception( "The source functionspace must have (halo >= 1) for pole treatment" ); } - if ( not matrix_free_ ) { - ATLAS_ASSERT( target_lonlat_ ); // TODO: implement setup with target_lonlat_fields_ as well (see execute_impl) - - idx_t inp_npts = source.size(); - idx_t out_npts = target_lonlat_.shape( 0 ); - - // return early if no output points on this partition reserve is called on - // the triplets but also during the sparseMatrix constructor. This won't - // work for empty matrices - if (out_npts == 0) { - return; - } + out_npts_ = 0; + if( target_lonlat_ ) { + convert_units_ = convert_units_multiplier(target_lonlat_); + out_npts_ = target_lonlat_.shape( 0 ); + } + else if( target_lonlat_fields_.empty() ) { + convert_units_ = convert_units_multiplier(target_lonlat_fields_[LON]); + out_npts_ = target_lonlat_fields_[0].shape(0); + } - auto lonlat = array::make_view( target_lonlat_ ); + if ( not matrix_free_ ) { + ATLAS_TRACE( "Precomputing interpolation matrix" ); - double convert_units = convert_units_multiplier( target_lonlat_ ); + std::vector failed_points; - auto triplets = kernel_->allocate_triplets( out_npts ); + auto triplets = kernel_->allocate_triplets( out_npts_ ); - const auto src_fs = functionspace::StructuredColumns( source ); - const auto src_dom = RectangularDomain( src_fs.grid().domain() ); - const double src_west = src_dom ? src_dom.xmin() : 0.; - const util::NormaliseLongitude normalise( src_west ); + using WorkSpace = typename Kernel::WorkSpace; + auto [west, east] = compute_src_west_east(source); + auto normalise = util::NormaliseLongitude(west); + auto interpolate_point = [&]( idx_t n, PointLonLat&& p, WorkSpace& workspace ) -> int { + try { + p.lon() = normalise(p.lon()); + kernel_->insert_triplets( n, p, triplets, workspace ); + return 0; + } + catch(const eckit::Exception& e) {} + if (verbose_) { + Log::error() << "Could not interpolate point " << n << " :\t" << p << std::endl; + } + return 1; + }; - ATLAS_TRACE_SCOPE( "Precomputing interpolation matrix" ) { - if ( target_ghost_ ) { - auto ghost = array::make_view( target_ghost_ ); - atlas_omp_parallel { - typename Kernel::WorkSpace workspace; - atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { - if ( not ghost( n ) ) { - PointLonLat p{normalise( lonlat( n, LON ) ) * convert_units, - lonlat( n, LAT ) * convert_units}; - kernel_->insert_triplets( n, p, triplets, workspace ); + auto interpolate_omp = [&]( idx_t out_npts, auto lonlat, auto ghost) { + atlas_omp_parallel { + WorkSpace workspace; + atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { + if( not ghost(n) ) { + if (interpolate_point(n, lonlat(n), workspace) != 0) { + atlas_omp_critical { + failed_points.emplace_back(n); + } } } } } - else { - atlas_omp_parallel { - typename Kernel::WorkSpace workspace; - atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { - PointLonLat p{normalise( lonlat( n, LON ) ) * convert_units, lonlat( n, LAT ) * convert_units}; - kernel_->insert_triplets( n, p, triplets, workspace ); - } + }; + + if ( target_lonlat_ ) { + auto lonlat_view = array::make_view( target_lonlat_ ); + auto lonlat = [lonlat_view, convert_units = convert_units_] (idx_t n) { + return PointLonLat{lonlat_view(n,LON) * convert_units, lonlat_view(n,LAT) * convert_units}; + }; + + if( out_npts_ != 0 ) { + if ( target_ghost_ ) { + auto ghost = array::make_view( target_ghost_ ); + interpolate_omp(out_npts_, lonlat, ghost); + } + else { + auto no_ghost = [](idx_t n) { return false; }; + interpolate_omp(out_npts_, lonlat, no_ghost); + } + } + handle_failed_points(*this, failed_points, lonlat); + } + else if ( not target_lonlat_fields_.empty() ) { + const auto lon = array::make_view( target_lonlat_fields_[LON] ); + const auto lat = array::make_view( target_lonlat_fields_[LAT] ); + auto lonlat = [lon, lat, convert_units = convert_units_] (idx_t n) { + return PointLonLat{lon(n) * convert_units, lat(n) * convert_units}; + }; + + if( out_npts_ != 0) { + if ( target_ghost_ ) { + auto ghost = array::make_view( target_ghost_ ); + interpolate_omp(out_npts_, lonlat, ghost); + } + else { + auto no_ghost = [](idx_t n) { return false; }; + interpolate_omp(out_npts_, lonlat, no_ghost); } } - // fill sparse matrix and return - Matrix A( out_npts, inp_npts, triplets ); + handle_failed_points(*this, failed_points, lonlat); + } + else { + ATLAS_NOTIMPLEMENTED; + } + + // fill sparse matrix + if( failed_points.empty() ) { + idx_t inp_npts = source.size(); + Matrix A( out_npts_, inp_npts, triplets ); setMatrix(A); } } @@ -256,12 +486,9 @@ template template void StructuredInterpolation2D::execute_impl( const Kernel& kernel, const FieldSet& src_fields, FieldSet& tgt_fields ) const { - const idx_t N = src_fields.size(); + using namespace structured2d; - const auto src_fs = functionspace::StructuredColumns( source() ); - const auto src_dom = RectangularDomain( src_fs.grid().domain() ); - const double src_west = src_dom ? src_dom.xmin() : 0.; - const util::NormaliseLongitude normalise( src_west ); + const idx_t N = src_fields.size(); std::vector > src_view; std::vector > tgt_view; @@ -272,65 +499,81 @@ void StructuredInterpolation2D::execute_impl( const Kernel& kernel, cons src_view.emplace_back( array::make_view( src_fields[i] ) ); tgt_view.emplace_back( array::make_view( tgt_fields[i] ) ); } - if ( target_lonlat_ ) { - double convert_units = convert_units_multiplier( target_lonlat_ ); - if ( target_ghost_ ) { - idx_t out_npts = target_lonlat_.shape( 0 ); - auto ghost = array::make_view( target_ghost_ ); - auto lonlat = array::make_view( target_lonlat_ ); + using WorkSpace = typename Kernel::WorkSpace; + + auto [west, east] = compute_src_west_east(source()); + util::NormaliseLongitude normalise{west}; + auto interpolate_point = [&]( idx_t n, PointLonLat&& p, WorkSpace& workspace ) -> int { + try { + p.lon() = normalise(p.lon()); + kernel.compute_stencil( p.lon(), p.lat(), workspace.stencil ); + kernel.compute_weights( p.lon(), p.lat(), workspace.stencil, workspace.weights ); + for ( idx_t i = 0; i < N; ++i ) { + kernel.interpolate( workspace.stencil, workspace.weights, src_view[i], tgt_view[i], n ); + } + return 0; + } + catch(const eckit::Exception& e) {} + if (verbose_) { + Log::error() << "Could not interpolate point " << n << " :\t" << p << std::endl; + } + return 1; + }; - atlas_omp_parallel { - typename Kernel::Stencil stencil; - typename Kernel::Weights weights; - atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { - if ( not ghost( n ) ) { - PointLonLat p{ normalise(lonlat( n, LON ) ) * convert_units, lonlat( n, LAT ) * convert_units}; - kernel.compute_stencil( p.lon(), p.lat(), stencil ); - kernel.compute_weights( p.lon(), p.lat(), stencil, weights ); - for ( idx_t i = 0; i < N; ++i ) { - kernel.interpolate( stencil, weights, src_view[i], tgt_view[i], n ); + std::vector failed_points; + + auto interpolate_omp = [&]( idx_t out_npts, auto lonlat, auto ghost) { + atlas_omp_parallel { + WorkSpace workspace; + atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { + if( not ghost(n) ) { + if (interpolate_point(n, lonlat(n), workspace) != 0) { + atlas_omp_critical { + failed_points.emplace_back(n); } } } } } - else { - idx_t out_npts = target_lonlat_.shape( 0 ); - const auto lonlat = array::make_view( target_lonlat_ ); + }; - atlas_omp_parallel { - typename Kernel::Stencil stencil; - typename Kernel::Weights weights; - atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { - PointLonLat p{ normalise(lonlat( n, LON )) * convert_units, lonlat( n, LAT ) * convert_units}; - kernel.compute_stencil( p.lon(), p.lat(), stencil ); - kernel.compute_weights( p.lon(), p.lat(), stencil, weights ); - for ( idx_t i = 0; i < N; ++i ) { - kernel.interpolate( stencil, weights, src_view[i], tgt_view[i], n ); - } - } + if ( target_lonlat_ ) { + auto lonlat_view = array::make_view( target_lonlat_ ); + auto lonlat = [lonlat_view, convert_units = convert_units_] (idx_t n) { + return PointLonLat{lonlat_view(n,LON) * convert_units, lonlat_view(n,LAT) * convert_units}; + }; + + if( out_npts_ != 0 ) { + if ( target_ghost_ ) { + auto ghost = array::make_view( target_ghost_ ); + interpolate_omp(out_npts_, lonlat, ghost); + } + else { + auto no_ghost = [](idx_t n) { return false; }; + interpolate_omp(out_npts_, lonlat, no_ghost); } } + handle_failed_points(*this, failed_points, lonlat); } else if ( not target_lonlat_fields_.empty() ) { - idx_t out_npts = target_lonlat_fields_[0].shape( 0 ); - const auto lon = array::make_view( target_lonlat_fields_[LON] ); - const auto lat = array::make_view( target_lonlat_fields_[LAT] ); - double convert_units = convert_units_multiplier( target_lonlat_fields_[LON] ); + const auto lon = array::make_view( target_lonlat_fields_[LON] ); + const auto lat = array::make_view( target_lonlat_fields_[LAT] ); + auto lonlat = [lon, lat, convert_units = convert_units_] (idx_t n) { + return PointLonLat{lon(n) * convert_units, lat(n) * convert_units}; + }; - atlas_omp_parallel { - typename Kernel::Stencil stencil; - typename Kernel::Weights weights; - atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { - PointLonLat p{normalise(lon( n )) * convert_units, lat( n ) * convert_units}; - kernel.compute_stencil( p.lon(), p.lat(), stencil ); - kernel.compute_weights( p.lon(), p.lat(), stencil, weights ); - for ( idx_t i = 0; i < N; ++i ) { - kernel.interpolate( stencil, weights, src_view[i], tgt_view[i], n ); - } + if( out_npts_ != 0) { + if ( target_ghost_ ) { + auto ghost = array::make_view( target_ghost_ ); + interpolate_omp(out_npts_, lonlat, ghost); + } + else { + auto no_ghost = [](idx_t n) { return false; }; + interpolate_omp(out_npts_, lonlat, no_ghost); } } + handle_failed_points(*this, failed_points, lonlat); } else { ATLAS_NOTIMPLEMENTED; diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc index 359749e1a..2c2b0a8d6 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc @@ -322,9 +322,9 @@ void StructuredInterpolation3D::execute_impl( const Kernel& kernel, cons typename Kernel::Weights weights; atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { for ( idx_t k = 0; k < out_nlev; ++k ) { - const double x = xcoords( n, k ) * convert_units; - const double y = ycoords( n, k ) * convert_units; - const double z = zcoords( n, k ); + double x = xcoords( n, k ) * convert_units; + double y = ycoords( n, k ) * convert_units; + double z = zcoords( n, k ); kernel.compute_stencil( x, y, z, stencil ); kernel.compute_weights( x, y, z, stencil, weights ); for ( idx_t i = 0; i < N; ++i ) { diff --git a/src/atlas/interpolation/method/structured/kernels/Cubic3DKernel.h b/src/atlas/interpolation/method/structured/kernels/Cubic3DKernel.h index ec6e1c088..55a4448b4 100644 --- a/src/atlas/interpolation/method/structured/kernels/Cubic3DKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/Cubic3DKernel.h @@ -73,13 +73,13 @@ class Cubic3DKernel { }; template - void compute_stencil(const double x, const double y, const double z, stencil_t& stencil) const { + void compute_stencil(double& x, const double y, const double z, stencil_t& stencil) const { horizontal_interpolation_.compute_stencil(x, y, stencil); vertical_interpolation_.compute_stencil(z, stencil); } template - void compute_weights(const double x, const double y, const double z, weights_t& weights) const { + void compute_weights(double x, double y, double z, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, z, stencil); compute_weights(x, y, z, stencil, weights); diff --git a/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h b/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h index 16fcbf1c6..da5235d6e 100644 --- a/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h @@ -72,12 +72,37 @@ class CubicHorizontalKernel { }; template - void compute_stencil(const double x, const double y, stencil_t& stencil) const { + void compute_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const { compute_horizontal_stencil_(x, y, stencil); + for (idx_t j = 0; j < stencil_width(); ++j) { + idx_t imin = stencil.i(0, j); + idx_t imax = stencil.i(stencil_width()-1, j); + if (imin < src_.i_begin_halo(stencil.j(j))) { + if (retry ) { + x += 360.; + compute_stencil(x, y, stencil, false); + } + else { + Log::error() << "Stencil out of bounds" << std::endl; + ATLAS_THROW_EXCEPTION("stencil out of bounds"); + } + } + if (imax >= src_.i_end_halo(stencil.j(j))) { + if (retry ) { + x -= 360.; + compute_stencil(x, y, stencil, false); + } + else { + Log::error() << "Stencil out of bounds" << std::endl; + ATLAS_THROW_EXCEPTION("Stencil out of bounds"); + } + } + } + } template - void compute_weights(const double x, const double y, weights_t& weights) const { + void compute_weights(double x, const double y, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); @@ -204,7 +229,7 @@ class CubicHorizontalKernel { } template - typename array_t::value_type operator()(const double x, const double y, const array_t& input) const { + typename array_t::value_type operator()(double x, double y, const array_t& input) const { Stencil stencil; compute_stencil(x, y, stencil); Weights weights; @@ -213,9 +238,9 @@ class CubicHorizontalKernel { } template - typename array_t::value_type interpolate(const PointLonLat& p, const array_t& input, WorkSpace& ws) const { - compute_stencil(p.lon(), p.lat(), ws.stencil); - compute_weights(p.lon(), p.lat(), ws.stencil, ws.weights); + typename array_t::value_type interpolate(PointXY p, const array_t& input, WorkSpace& ws) const { + compute_stencil(p.x(), p.y(), ws.stencil); + compute_weights(p.x(), p.y(), ws.stencil, ws.weights); return interpolate(ws.stencil, ws.weights, input); } @@ -239,8 +264,8 @@ class CubicHorizontalKernel { insert_triplets(row, p.x(), p.y(), triplets, ws); } - void insert_triplets(const idx_t row, const double x, const double y, Triplets& triplets, WorkSpace& ws) const { - compute_horizontal_stencil_(x, y, ws.stencil); + void insert_triplets(const idx_t row, double x, double y, Triplets& triplets, WorkSpace& ws) const { + compute_stencil(x, y, ws.stencil); compute_weights(x, y, ws.stencil, ws.weights); const auto& wj = ws.weights.weights_j; diff --git a/src/atlas/interpolation/method/structured/kernels/Linear3DKernel.h b/src/atlas/interpolation/method/structured/kernels/Linear3DKernel.h index e1d9ad99b..5efee6496 100644 --- a/src/atlas/interpolation/method/structured/kernels/Linear3DKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/Linear3DKernel.h @@ -65,13 +65,13 @@ class Linear3DKernel { }; template - void compute_stencil(const double x, const double y, const double z, stencil_t& stencil) const { + void compute_stencil(double& x, const double y, const double z, stencil_t& stencil) const { horizontal_interpolation_.compute_stencil(x, y, stencil); vertical_interpolation_.compute_stencil(z, stencil); } template - void compute_weights(const double x, const double y, const double z, weights_t& weights) const { + void compute_weights(double x, double y, const double z, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, z, stencil); compute_weights(x, y, z, stencil, weights); diff --git a/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h b/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h index 858606272..cccab6d0c 100644 --- a/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h @@ -64,12 +64,37 @@ class LinearHorizontalKernel { }; template - void compute_stencil(const double x, const double y, stencil_t& stencil) const { + void compute_stencil(double& x, double y, stencil_t& stencil, bool retry = true) const { compute_horizontal_stencil_(x, y, stencil); + for (idx_t j = 0; j < stencil_width(); ++j) { + idx_t imin = stencil.i(0, j); + idx_t imax = stencil.i(stencil_width()-1, j); + if (imin < src_.i_begin_halo(stencil.j(j))) { + if (retry ) { + x += 360.; + compute_stencil(x, y, stencil, false); + } + else { + Log::error() << "Stencil out of bounds" << std::endl; + ATLAS_THROW_EXCEPTION("stencil out of bounds"); + } + } + if (imax >= src_.i_end_halo(stencil.j(j))) { + if (retry ) { + x -= 360.; + compute_stencil(x, y, stencil, false); + } + else { + Log::error() << "Stencil out of bounds" << std::endl; + ATLAS_THROW_EXCEPTION("Stencil out of bounds"); + } + } + } } + template - void compute_weights(const double x, const double y, weights_t& weights) const { + void compute_weights(double x, double y, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); @@ -156,18 +181,18 @@ class LinearHorizontalKernel { } template - typename array_t::value_type operator()(const double x, const double y, const array_t& input) const { + typename array_t::value_type operator()(double x, double y, const array_t& input) const { Stencil stencil; - compute_stencil(x, y, stencil); Weights weights; + compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); return interpolate(stencil, weights, input); } template - typename array_t::value_type interpolate(const PointLonLat& p, const array_t& input, WorkSpace& ws) const { - compute_stencil(p.lon(), p.lat(), ws.stencil); - compute_weights(p.lon(), p.lat(), ws.stencil, ws.weights); + typename array_t::value_type interpolate(PointXY p, const array_t& input, WorkSpace& ws) const { + compute_stencil(p.x(), p.y(), ws.stencil); + compute_weights(p.x(), p.y(), ws.stencil, ws.weights); return interpolate(ws.stencil, ws.weights, input); } @@ -191,8 +216,8 @@ class LinearHorizontalKernel { insert_triplets(row, p.x(), p.y(), triplets, ws); } - void insert_triplets(const idx_t row, const double x, const double y, Triplets& triplets, WorkSpace& ws) const { - compute_horizontal_stencil_(x, y, ws.stencil); + void insert_triplets(const idx_t row, double x, double y, Triplets& triplets, WorkSpace& ws) const { + compute_stencil(x, y, ws.stencil, true); compute_weights(x, y, ws.stencil, ws.weights); const auto& wj = ws.weights.weights_j; diff --git a/src/atlas/interpolation/method/structured/kernels/QuasiCubic3DKernel.h b/src/atlas/interpolation/method/structured/kernels/QuasiCubic3DKernel.h index 343fc6d19..aeb7c3424 100644 --- a/src/atlas/interpolation/method/structured/kernels/QuasiCubic3DKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/QuasiCubic3DKernel.h @@ -89,13 +89,13 @@ class QuasiCubic3DKernel { }; template - void compute_stencil(const double x, const double y, const double z, stencil_t& stencil) const { + void compute_stencil(double& x, const double y, const double z, stencil_t& stencil) const { quasi_cubic_horizontal_interpolation_.compute_stencil(x, y, stencil); vertical_interpolation_.compute_stencil(z, stencil); } template - void compute_weights(const double x, const double y, const double z, weights_t& weights) const { + void compute_weights(double x, double y, double z, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, z, stencil); compute_weights(x, y, z, stencil, weights); diff --git a/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h b/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h index 4a3d8bd75..e9280bd76 100644 --- a/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h @@ -74,12 +74,37 @@ class QuasiCubicHorizontalKernel { }; template - void compute_stencil(const double x, const double y, stencil_t& stencil) const { + void compute_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const { compute_horizontal_stencil_(x, y, stencil); + for (idx_t j = 0; j < stencil_width(); ++j) { + idx_t imin = stencil.i(0, j); + idx_t imax = stencil.i(stencil_width()-1, j); + if (imin < src_.i_begin_halo(stencil.j(j))) { + if (retry ) { + x += 360.; + compute_stencil(x, y, stencil, false); + } + else { + Log::error() << "Stencil out of bounds" << std::endl; + ATLAS_THROW_EXCEPTION("stencil out of bounds"); + } + } + if (imax >= src_.i_end_halo(stencil.j(j))) { + if (retry ) { + x -= 360.; + compute_stencil(x, y, stencil, false); + } + else { + Log::error() << "Stencil out of bounds" << std::endl; + ATLAS_THROW_EXCEPTION("Stencil out of bounds"); + } + } + } + } template - void compute_weights(const double x, const double y, weights_t& weights) const { + void compute_weights(double x, double y, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); @@ -255,7 +280,7 @@ class QuasiCubicHorizontalKernel { } template - typename array_t::value_type operator()(const double x, const double y, const array_t& input) const { + typename array_t::value_type operator()(double x, double y, const array_t& input) const { Stencil stencil; compute_stencil(x, y, stencil); Weights weights; @@ -264,9 +289,9 @@ class QuasiCubicHorizontalKernel { } template - typename array_t::value_type interpolate(const PointLonLat& p, const array_t& input, WorkSpace& ws) const { - compute_stencil(p.lon(), p.lat(), ws.stencil); - compute_weights(p.lon(), p.lat(), ws.stencil, ws.weights); + typename array_t::value_type interpolate(PointXY p, const array_t& input, WorkSpace& ws) const { + compute_stencil(p.x(), p.y(), ws.stencil); + compute_weights(p.x(), p.y(), ws.stencil, ws.weights); return interpolate(ws.stencil, ws.weights, input); } @@ -290,8 +315,8 @@ class QuasiCubicHorizontalKernel { insert_triplets(row, p.x(), p.y(), triplets, ws); } - void insert_triplets(const idx_t row, const double x, const double y, Triplets& triplets, WorkSpace& ws) const { - compute_horizontal_stencil_(x, y, ws.stencil); + void insert_triplets(const idx_t row, double x, double y, Triplets& triplets, WorkSpace& ws) const { + compute_stencil(x, y, ws.stencil); compute_weights(x, y, ws.stencil, ws.weights); const auto& wj = ws.weights.weights_j; diff --git a/src/tests/interpolation/test_interpolation_structured2D.cc b/src/tests/interpolation/test_interpolation_structured2D.cc index c5847f03d..4a668dd5d 100644 --- a/src/tests/interpolation/test_interpolation_structured2D.cc +++ b/src/tests/interpolation/test_interpolation_structured2D.cc @@ -51,6 +51,7 @@ static Config scheme() { scheme.set("halo", 2); } scheme.set("name", scheme_str); + scheme.set("verbose",eckit::Resource("--verbose",false)); return scheme; } From edecc726daea74d0d8bbf21558c1cf10dec62ac7 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 28 Mar 2023 10:57:27 +0000 Subject: [PATCH 16/19] Fix compilation with nvidia/22.1 --- .../structured/StructuredInterpolation2D.tcc | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index 8f9ae1392..e18cd06bd 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -41,6 +41,22 @@ #include "atlas/util/NormaliseLongitude.h" #include "atlas/util/Point.h" +// With nvidia/22.1, following error is encountered during interpolation. +// NVC++-S-0000-Internal compiler error. BAD sptr in var_refsym 0 +// It is observed to have been fixed with nvidia/22.11 +// Disabling OpenMP in this routine for nvidia < 22.11 seems to fix things +#if defined(__NVCOMPILER) +#if (__NVCOMPILER_MAJOR__*100) + __NVCOMPILER_MINOR__ < 2211 +#warning "Disabled OpenMP for StructuredInterpolation2D due to internal compiler error" +#undef atlas_omp_parallel +#define atlas_omp_parallel +#undef atlas_omp_for +#define atlas_omp_for for +#endif +#endif + + + namespace atlas { namespace interpolation { namespace method { @@ -365,7 +381,7 @@ void StructuredInterpolation2D::setup( const FunctionSpace& source ) { return 1; }; - auto interpolate_omp = [&]( idx_t out_npts, auto lonlat, auto ghost) { + auto interpolate_omp = [&failed_points,interpolate_point]( idx_t out_npts, auto lonlat, auto ghost) { atlas_omp_parallel { WorkSpace workspace; atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { @@ -523,7 +539,7 @@ void StructuredInterpolation2D::execute_impl( const Kernel& kernel, cons std::vector failed_points; - auto interpolate_omp = [&]( idx_t out_npts, auto lonlat, auto ghost) { + auto interpolate_omp = [&failed_points,interpolate_point]( idx_t out_npts, auto lonlat, auto ghost) { atlas_omp_parallel { WorkSpace workspace; atlas_omp_for( idx_t n = 0; n < out_npts; ++n ) { From 62d08a024341799104486e032732153c4280d4ab Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 28 Mar 2023 09:10:36 +0000 Subject: [PATCH 17/19] Avoid some compiler warnings with nvhpc --- atlas_io/src/atlas_io/types/scalar.cc | 2 +- src/atlas/functionspace/detail/BlockStructuredColumns.h | 2 ++ .../interpolation/method/cubedsphere/CubedSphereBilinear.h | 5 +++-- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/atlas_io/src/atlas_io/types/scalar.cc b/atlas_io/src/atlas_io/types/scalar.cc index 61448c82c..c4ce6a80a 100644 --- a/atlas_io/src/atlas_io/types/scalar.cc +++ b/atlas_io/src/atlas_io/types/scalar.cc @@ -14,7 +14,7 @@ #endif // GNU C++ compiler (version 11) should not try to optimize this code -#ifdef __GNUC__ +#if defined(__GNUC__) && !defined(__NVCOMPILER) #pragma GCC optimize("O0") #endif diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.h b/src/atlas/functionspace/detail/BlockStructuredColumns.h index c34a0aaa7..6e4919fe3 100644 --- a/src/atlas/functionspace/detail/BlockStructuredColumns.h +++ b/src/atlas/functionspace/detail/BlockStructuredColumns.h @@ -83,8 +83,10 @@ class BlockStructuredColumns : public FunctionSpaceImpl { Field createField(const eckit::Configuration&) const override; Field createField(const Field&, const eckit::Configuration&) const override; + using FunctionSpaceImpl::scatter; void scatter(const FieldSet&, FieldSet&) const override; void scatter(const Field&, Field&) const override; + using FunctionSpaceImpl::gather; void gather(const FieldSet&, FieldSet&) const override; void gather(const Field&, Field&) const override; diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h index 166ec8c61..ed7dce143 100644 --- a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h +++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h @@ -39,8 +39,9 @@ class CubedSphereBilinear : public Method { virtual const FunctionSpace& target() const override { return target_; } private: - virtual void do_setup(const FunctionSpace& source, const FunctionSpace& target) override; - virtual void do_setup(const Grid& source, const Grid& target, const Cache&) override; + using Method::do_setup; + void do_setup(const FunctionSpace& source, const FunctionSpace& target) override; + void do_setup(const Grid& source, const Grid& target, const Cache&) override; FunctionSpace source_; FunctionSpace target_; From 24949b608f8699f661664ebd083dc09aaec379af Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Tue, 7 Mar 2023 16:52:59 +0100 Subject: [PATCH 18/19] Add MDPI functions from Valcke et al (2022) paper for gauging the quality of Atlas interpolations --- src/atlas/CMakeLists.txt | 2 + src/atlas/util/function/MDPI_functions.cc | 123 ++++++++++++++++++ src/atlas/util/function/MDPI_functions.h | 37 ++++++ .../atlas-conservative-interpolation.cc | 17 +++ 4 files changed, 179 insertions(+) create mode 100644 src/atlas/util/function/MDPI_functions.cc create mode 100644 src/atlas/util/function/MDPI_functions.h diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 938ba1ab2..3ca7f0b55 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -757,6 +757,8 @@ util/detail/BlackMagic.h util/detail/Cache.h util/detail/Debug.h util/detail/KDTree.h +util/function/MDPI_functions.h +util/function/MDPI_functions.cc util/function/SolidBodyRotation.h util/function/SolidBodyRotation.cc util/function/SphericalHarmonic.h diff --git a/src/atlas/util/function/MDPI_functions.cc b/src/atlas/util/function/MDPI_functions.cc new file mode 100644 index 000000000..121d2602d --- /dev/null +++ b/src/atlas/util/function/MDPI_functions.cc @@ -0,0 +1,123 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include + +#include "atlas/util/Constants.h" +#include "atlas/util/Earth.h" + +#include "atlas/util/function/MDPI_functions.h" + +namespace atlas { + +namespace util { + +namespace function { + +double MDPI_sinusoid(double lon, double lat) { + lon *= Constants::degreesToRadians(); + lat *= Constants::degreesToRadians(); + constexpr double length = 1.2 * M_PI; + + return 2. - std::cos(M_PI * std::acos(std::cos(lon) * std::cos(lat)) / length); +} + +double MDPI_harmonic(double lon, double lat) { + lon *= Constants::degreesToRadians(); + lat *= Constants::degreesToRadians(); + + return 2. + std::pow(std::sin(2. * lat), 16) * std::cos(16. * lon); +} + +double MDPI_vortex(double lon, double lat) { + lon *= Constants::degreesToRadians(); + lat *= Constants::degreesToRadians(); + constexpr double dLon0 = 5.5; + constexpr double dLat0 = 0.2; + constexpr double dR0 = 3.0; + constexpr double dD = 5.0; + constexpr double dT = 6.0; + constexpr double length = 1.2 * M_PI; + + auto sqr = [](const double x) { return x * x; }; + auto sech = [](const double x) { return 1. / std::cosh(x); }; + + double dSinC = std::sin(dLat0); + double dCosC = std::cos(dLat0); + double dCosT = std::cos(lat); + double dSinT = std::sin(lat); + double dTrm = dCosT * std::cos(lon - dLon0); + double dX = dSinC * dTrm - dCosC * dSinT; + double dY = dCosT * std::sin(lon - dLon0); + double dZ = dSinC * dSinT + dCosC * dTrm; + + double dlon = std::atan2(dY, dX); + double dlat = std::asin(dZ); + + double dRho = dR0 * std::cos(dlat); + double dVt = 1.5 * std::sqrt(3.) * sqr(sech(dRho)) * std::tanh(dRho); + double dOmega; + if (dRho == 0.) { + dOmega = 0.; + } + else { + dOmega = dVt / dRho; + } + + return 2. * (1. + std::tanh(dRho / dD * std::sin(dlon - dOmega * dT))); +} + +double MDPI_gulfstream(double lon, double lat) { + constexpr double d2r = Constants::degreesToRadians(); + + auto sqr = [](const double x) { return x * x; }; + auto sech = [](const double x) { return 1. / std::cosh(x); }; + + constexpr double length = 1.2 * M_PI; + constexpr double gf_coef = 1.0; // Coefficient for a Gult Stream term (0.0 = no Gulf Stream) + constexpr double gf_ori_lon = -80.0 * d2r; // Origin of the Gulf Stream (longitude in deg) + constexpr double gf_ori_lat = 25.0 * d2r; // Origin of the Gulf Stream (latitude in deg) + constexpr double gf_end_lon = -1.8 * d2r; // End of the Gulf Stream (longitude in deg) + constexpr double gf_end_lat = 50.0 * d2r; // End of the Gulf Stream (latitude in deg) + constexpr double gf_dmp_lon = -25.5 * d2r; // Point of the Gulf Stream decrease (longitude in deg) + constexpr double gf_dmp_lat = -55.5 * d2r; // Point of the Gulf Stream decrease (latitude in deg) + + double dr0 = std::sqrt(sqr(gf_end_lon - gf_ori_lon) + sqr(gf_end_lat - gf_ori_lat)); + double dr1 = std::sqrt(sqr(gf_dmp_lon - gf_ori_lon) + sqr(gf_dmp_lat - gf_ori_lat)); + + double gf_per_lon = [lon,d2r]() { + double gf_per_lon = lon - 180.; + while (gf_per_lon > 180.) { + gf_per_lon -= 360.; + } + while (gf_per_lon < -180.) { + gf_per_lon += 360.; + } + return gf_per_lon * d2r; + }(); + double dx = gf_per_lon - gf_ori_lon; + double dy = lat * d2r - gf_ori_lat; + double dr = std::sqrt(sqr(dx) + sqr(dy)); + double dth = std::atan2(dy, dx); + double dc = 1.3 * gf_coef; + if (dr > dr0) { + dc = 0.; + } + if (dr > dr1) { + dc *= std::cos(M_PI_2 * (dr - dr1)/(dr0 - dr1)); + } + + double background_func = MDPI_sinusoid(lon, lat); + return background_func + dc * (std::max(1000. * std::sin(0.4 * (0.5 * dr + dth) + 0.007 * std::cos(50. * dth) + 0.37 * M_PI), 999.) - 999.); +} + +} // namespace function +} // namespace util +} // namespace atlas diff --git a/src/atlas/util/function/MDPI_functions.h b/src/atlas/util/function/MDPI_functions.h new file mode 100644 index 000000000..b60544536 --- /dev/null +++ b/src/atlas/util/function/MDPI_functions.h @@ -0,0 +1,37 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation nor + * does it submit to any jurisdiction. + */ + + +#pragma once + +namespace atlas { + +namespace util { + +namespace function { + +/// \brief Analytic functions from MDPI paper on regridding +/// +/// \detailed The formula is found in +/// "Benchmarking Regridding Libraries Used in Earth System Modelling" +/// by Sophie Valcke, Andreas Piacentini, Gabriel Jonville, MDPI 2022 +/// as the sinusoid analytical function in Sec 2.1.2. +/// The longitude (lon) and latitude (lat) are assumed to be in degrees, +/// +double MDPI_sinusoid(double lon, double lat); +double MDPI_harmonic(double lon, double lat); +double MDPI_vortex(double lon, double lat); +double MDPI_gulfstream(double lon, double lat); + +} // namespace function + +} // namespace util + +} // namespace atlas diff --git a/src/sandbox/interpolation/atlas-conservative-interpolation.cc b/src/sandbox/interpolation/atlas-conservative-interpolation.cc index afe570c92..48b6e71eb 100644 --- a/src/sandbox/interpolation/atlas-conservative-interpolation.cc +++ b/src/sandbox/interpolation/atlas-conservative-interpolation.cc @@ -45,6 +45,7 @@ #include "atlas/output/Gmsh.h" #include "atlas/runtime/AtlasTool.h" #include "atlas/util/Config.h" +#include "atlas/util/function/MDPI_functions.h" #include "atlas/util/function/SolidBodyRotation.h" #include "atlas/util/function/SphericalHarmonic.h" #include "atlas/util/function/VortexRollup.h" @@ -155,6 +156,22 @@ std::function get_init(const AtlasTool::Args& args) util::function::SolidBodyRotation sbr(beta); return [sbr](const PointLonLat& p) { return sbr.windMagnitude(p.lon(), p.lat()); }; } + else if (init == "MDPI_sinusoid") { + auto sbr = util::function::MDPI_sinusoid; + return [sbr](const PointLonLat& p) { return sbr(p.lon(), p.lat()); }; + } + else if (init == "MDPI_harmonic") { + auto sbr = util::function::MDPI_harmonic; + return [sbr](const PointLonLat& p) { return sbr(p.lon(), p.lat()); }; + } + else if (init == "MDPI_vortex") { + auto sbr = util::function::MDPI_vortex; + return [sbr](const PointLonLat& p) { return sbr(p.lon(), p.lat()); }; + } + else if (init == "MDPI_gulfstream") { + auto sbr = util::function::MDPI_gulfstream; + return [sbr](const PointLonLat& p) { return sbr(p.lon(), p.lat()); }; + } else { if (args.has("init")) { Log::error() << "Bad value for \"init\": \"" << init << "\" not recognised." << std::endl; From c85c9f401e8d2daa2c7748319d28236386b6b034 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 30 Mar 2023 16:30:03 +0200 Subject: [PATCH 19/19] Version 0.33.0 --- CHANGELOG.md | 12 ++++++++++++ VERSION | 2 +- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 09676f01e..e04f14468 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,17 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] +## [0.33.0] - 2023-04-03 +### Added +- Add support for StructuredPartitionPolygon with halo > 0 +- Add information on atlas having PROJ support + +### Changed +- C++17 standard is now a requirement + +### Fixed +Fix StructuredInterpolation2D with retry for failed stencils + ## [0.32.1] - 2023-02-09 ### Added - Added (lon, lat) to (alpha, beta) transforms to cubed sphere projection @@ -432,6 +443,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## 0.13.0 - 2018-02-16 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop +[0.33.0]: https://github.com/ecmwf/atlas/compare/0.32.1...0.33.0 [0.32.1]: https://github.com/ecmwf/atlas/compare/0.32.0...0.32.1 [0.32.0]: https://github.com/ecmwf/atlas/compare/0.31.1...0.32.0 [0.31.1]: https://github.com/ecmwf/atlas/compare/0.31.0...0.31.1 diff --git a/VERSION b/VERSION index fd9620c08..be386c9ed 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.32.1 +0.33.0