Merge branch 'develop' into feature/interpolation_comparisons

.github/ci-config.yml

@@ -1,5 +1,6 @@
 dependencies: |
   ecmwf/ecbuild
   ecmwf/eckit
+  ecmwf/fckit
 dependency_branch: develop
 parallelism_factor: 8

.github/ci-hpc-config.yml

@@ -12,6 +12,7 @@ mpi_on:
     dependencies:
       - ecmwf/ecbuild@develop
       - ecmwf/eckit@develop
+      - ecmwf/fckit@develop
     parallel: 64
     ntasks: 16
 
@@ -24,4 +25,5 @@ mpi_off:
     dependencies:
       - ecmwf/ecbuild@develop
       - ecmwf/eckit@develop
+      - ecmwf/fckit@develop
     parallel: 64

.github/workflows/build.yml

@@ -194,6 +194,12 @@ jobs:
         ${ATLAS_TOOLS}/install-fftw.sh --version 3.3.10 --prefix ${DEPS_DIR}/fftw
         echo "FFTW_ROOT=${DEPS_DIR}/fftw" >> $GITHUB_ENV
 
+    - name: Install Qhull
+      shell: bash -eux {0}
+      run: |
+        ${ATLAS_TOOLS}/install-qhull.sh --version 8.1-alpha3 --prefix ${DEPS_DIR}/qhull
+        echo "Qhull_ROOT=${DEPS_DIR}/qhull" >> $GITHUB_ENV
+
     - name: Install LZ4
       if: "!contains( matrix.compiler, 'nvhpc' )"
       run: |

cmake/features/ECTRANS.cmake

@@ -2,13 +2,16 @@ if( atlas_HAVE_ATLAS_TRANS )
 
 ### trans ...
 
+if( NOT DEFINED ATLAS_ENABLE_ECTRANS AND DEFINED ATLAS_ENABLE_TRANS )
+  ecbuild_warn("Atlas option ATLAS_ENABLE_TRANS is deprecated in favour of ATLAS_ENABLE_ECTRANS")
+  set( ATLAS_ENABLE_ECTRANS ${ATLAS_ENABLE_TRANS} )
+endif()
 if( NOT DEFINED ENABLE_ECTRANS AND DEFINED ENABLE_TRANS )
   ecbuild_warn("Atlas option ENABLE_TRANS is deprecated in favour of ENABLE_ECTRANS")
   set( ENABLE_ECTRANS ${ENABLE_TRANS} )
 endif()
-if( NOT DEFINED ATLAS_ENABLE_ECTRANS AND DEFINED ATLAS_ENABLE_TRANS )
-  ecbuild_warn("Atlas option ATLAS_ENABLE_TRANS is deprecated in favour of ATLAS_ENABLE_ECTRANS")
-  set( ATLAS_ENABLE_TRANS ${ENABLE_TRANS} )
+if( DEFINED ATLAS_ENABLE_ECTRANS )
+  set( ENABLE_ECTRANS ${ATLAS_ENABLE_ECTRANS} )
 endif()
 
 set( atlas_HAVE_PACKAGE_ECTRANS 0 )
@@ -32,4 +35,4 @@ ecbuild_add_option( FEATURE ECTRANS
                     DESCRIPTION "Support for IFS spectral transforms"
                     CONDITION atlas_HAVE_ATLAS_FUNCTIONSPACE AND transi_FOUND )
 
-endif()
+endif()

cmake/features/TESSELATION.cmake

@@ -1,16 +1,32 @@
 if( atlas_HAVE_ATLAS_FUNCTIONSPACE )
 ### tesselation ...
 
-set(Boost_USE_MULTITHREADED      ON )
-
 ecbuild_add_option( FEATURE TESSELATION
+                    DESCRIPTION "Support for unstructured mesh generation"
+                    CONDITION atlas_HAVE_ATLAS_FUNCTIONSPACE
+                    REQUIRED_PACKAGES "Qhull" )
+if(HAVE_TESSELATION)
+  set(QHULL_LIBRARIES Qhull::qhullcpp Qhull::qhull_r)
+  set(atlas_HAVE_QHULL 1)
+else()
+  set(atlas_HAVE_QHULL 0)
+endif()
+
+### NOTE
+#
+# CGAL is deprecated as TESSELATION backend. Qhull is to be used instead.
+# To use CGAL regardless, turn ON CGAL feature (-DENABLE_CGAL=ON)
+
+set(Boost_USE_MULTITHREADED ON )
+ecbuild_add_option( FEATURE CGAL
+                    DEFAULT OFF
                     DESCRIPTION "Support for unstructured mesh generation"
                     CONDITION atlas_HAVE_ATLAS_FUNCTIONSPACE
                     REQUIRED_PACKAGES
-                      "CGAL QUIET"
+                      "CGAL"
                       "Boost VERSION 1.45.0 QUIET" )
 
-if( HAVE_TESSELATION )
+if( HAVE_CGAL )
     list( APPEND CGAL_INCLUDE_DIRS ${Boost_INCLUDE_DIRS} )
     if ( TARGET CGAL::CGAL )
       list( APPEND CGAL_LIBRARIES CGAL::CGAL ${CGAL_3RD_PARTY_LIBRARIES} ${GMP_LIBRARIES} ${MPFR_LIBRARIES} ${Boost_THREAD_LIBRARY} ${Boost_SYSTEM_LIBRARY} )
@@ -27,8 +43,8 @@ if( HAVE_TESSELATION )
     endif()
 endif()
 
-if( NOT HAVE_TESSELATION )
+if( NOT HAVE_CGAL )
     unset( CGAL_LIBRARIES )
     unset( CGAL_INCLUDE_DIRS )
 endif()
-endif()
+endif()

src/CMakeLists.txt

@@ -24,6 +24,18 @@ else()
   set( atlas_HAVE_TESSELATION 0 )
 endif()
 
+if( atlas_HAVE_QHULL )
+  set( atlas_HAVE_QHULL 1 )
+else()
+  set( atlas_HAVE_QHULL 0 )
+endif()
+
+if( atlas_HAVE_CGAL )
+  set( atlas_HAVE_CGAL 1 )
+else()
+  set( atlas_HAVE_CGAL 0 )
+endif()
+
 if( atlas_HAVE_PROJ )
   set( atlas_HAVE_PROJ 1 )
 else()

src/atlas/CMakeLists.txt

@@ -341,6 +341,10 @@ mesh/Connectivity.cc
 mesh/Connectivity.h
 mesh/ElementType.cc
 mesh/ElementType.h
+mesh/elementtypes/Classification.h
+mesh/elementtypes/Triangle.h
+mesh/elementtypes/Quadrilateral.h
+mesh/elementtypes/Pentagon.h
 mesh/Elements.cc
 mesh/Elements.h
 mesh/Halo.cc
@@ -815,6 +819,10 @@ util/vector.h
 util/mdspan.h
 util/detail/mdspan/mdspan.hpp
 util/VectorOfAbstract.h
+util/QhullSphericalTriangulation.h
+util/QhullSphericalTriangulation.cc
+util/CGALSphericalTriangulation.h
+util/CGALSphericalTriangulation.cc
 util/detail/Cache.h
 util/detail/KDTree.h
 util/function/MDPI_functions.h
@@ -928,6 +936,7 @@ ecbuild_add_library( TARGET atlas
     ${CGAL_LIBRARIES}
     ${FFTW_LIBRARIES}
     ${PROJ_LIBRARIES}
+    ${QHULL_LIBRARIES}
 
   PUBLIC_LIBS
     eckit

src/atlas/functionspace/PointCloud.cc

@@ -70,9 +70,22 @@ PointCloud::PointCloud(const Field& lonlat, const Field& ghost): lonlat_(lonlat)
     setupHaloExchange();
 }
 
-PointCloud::PointCloud(const FieldSet & flds): lonlat_(flds["lonlat"]),
-  ghost_(flds["ghost"]), remote_index_(flds["remote_index"]), partition_(flds["partition"]) {
-    setupHaloExchange();
+PointCloud::PointCloud(const FieldSet & flds): lonlat_(flds["lonlat"]) {
+    if (flds.has("ghost")) {
+        ghost_ = flds["ghost"];
+    }
+    if (flds.has("remote_index")) {
+        remote_index_ = flds["remote_index"];
+    }
+    if (flds.has("partition")) {
+        partition_ = flds["partition"];
+    }
+    if (flds.has("global_index")) {
+        global_index_ = flds["global_index"];
+    }
+    if( ghost_ && remote_index_ && partition_ ) {
+        setupHaloExchange();
+    }
 }
 
 PointCloud::PointCloud(const Grid& grid) {

src/atlas/functionspace/PointCloud.h

@@ -56,6 +56,7 @@ class PointCloud : public functionspace::FunctionSpaceImpl {
     const Field& vertical() const { return vertical_; }
     Field ghost() const override;
     Field remote_index() const override { return remote_index_; }
+    Field global_index() const override { return global_index_; }
     virtual idx_t size() const override { return lonlat_.shape(0); }
 
     using FunctionSpaceImpl::createField;
@@ -154,6 +155,7 @@ class PointCloud : public functionspace::FunctionSpaceImpl {
     Field vertical_;
     mutable Field ghost_;
     Field remote_index_;
+    Field global_index_;
     Field partition_;
     std::unique_ptr<parallel::HaloExchange> halo_exchange_;
     idx_t levels_{0};

src/atlas/interpolation/Cache.h

@@ -108,6 +108,7 @@ class MatrixCache final : public Cache {
 public:
     MatrixCache() = default;
     MatrixCache(const Cache& c);
+    MatrixCache(const MatrixCache& c) : MatrixCache(Cache(c)) {}
     MatrixCache(Matrix&& m);
     MatrixCache(std::shared_ptr<const Matrix> m, const std::string& uid = "");
     MatrixCache(const Matrix* m);
@@ -146,6 +147,7 @@ class IndexKDTreeCache final : public Cache {
 public:
     IndexKDTreeCache() = default;
     IndexKDTreeCache(const Cache& c);
+    IndexKDTreeCache(const IndexKDTreeCache& c) : IndexKDTreeCache(Cache(c)) {}
     IndexKDTreeCache(const IndexKDTree&);
     IndexKDTreeCache(const Interpolation&);
     operator bool() const;

src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc

@@ -105,18 +105,6 @@ namespace {
         return out.str();
     };
 
-    inline auto compute_src_west_east(functionspace::StructuredColumns src_fs) {
-        const auto& src_polygon = src_fs.polygon(0);
-        double src_west = std::numeric_limits<double>::max();
-        double src_east = std::numeric_limits<double>::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));
@@ -182,13 +170,11 @@ namespace {
             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<mpi::size(); ++p ) {
                 if( p == my_rank && failed_points.size() ) {
                     Log::error() << "Failed to interpolate " << failed_points.size() << " points on rank " << p
-                                 << " after retrying with normalisation to (west=" << src_west-360 << ") and (east=" << src_west+360 << "). "
-                                 << "See " << filename("atlas_target_failed_points_p%p.json") ;
+                                 << ".\nSee " << filename("atlas_target_failed_points_p%p.json") ;
                     if( failed_points.size() <= 20 ) {
                         Log::error() << " :\n" << to_str(failed_points, lonlat);
                     }
@@ -366,11 +352,8 @@ void StructuredInterpolation2D<Kernel>::setup( const FunctionSpace& source ) {
         auto triplets = kernel_->allocate_triplets( out_npts_ );
 
         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;
             }
@@ -518,13 +501,11 @@ void StructuredInterpolation2D<Kernel>::execute_impl( const Kernel& kernel, cons
 
     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 );
+            kernel.make_valid_stencil( p.lon(), p.lat(), workspace.stencil );
             for ( idx_t i = 0; i < N; ++i ) {
                 kernel.interpolate( workspace.stencil, workspace.weights, src_view[i], tgt_view[i], n );
             }

src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h

@@ -72,15 +72,20 @@ class CubicHorizontalKernel {
     };
 
     template <typename stencil_t>
-    void compute_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const {
+    void compute_stencil(const double x, const double y, stencil_t& stencil) const {
         compute_horizontal_stencil_(x, y, stencil);
+    }
+
+    template <typename stencil_t>
+    void make_valid_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const {
         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);
+                    compute_stencil(x, y, stencil);
+                    return make_valid_stencil(x, y, stencil, false);
                 }
                 else {
                     Log::error() << "Stencil out of bounds" << std::endl;
@@ -90,7 +95,8 @@ class CubicHorizontalKernel {
             if (imax >= src_.i_end_halo(stencil.j(j))) {
                 if (retry ) {
                     x -= 360.;
-                    compute_stencil(x, y, stencil, false);
+                    compute_stencil(x, y, stencil);
+                    return make_valid_stencil(x, y, stencil, false);
                 }
                 else {
                     Log::error() << "Stencil out of bounds" << std::endl;
@@ -102,7 +108,7 @@ class CubicHorizontalKernel {
     }
 
     template <typename weights_t>
-    void compute_weights(double x, const double y, weights_t& weights) const {
+    void compute_weights(const double x, const double y, weights_t& weights) const {
         Stencil stencil;
         compute_stencil(x, y, stencil);
         compute_weights(x, y, stencil, weights);
@@ -234,13 +240,15 @@ class CubicHorizontalKernel {
         compute_stencil(x, y, stencil);
         Weights weights;
         compute_weights(x, y, stencil, weights);
+        make_valid_stencil(x, y, stencil);
         return interpolate(stencil, weights, input);
     }
 
     template <typename array_t>
     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);
+        make_valid_stencil(p.x(), p.y(), ws.stencil);
         return interpolate(ws.stencil, ws.weights, input);
     }
 
@@ -267,6 +275,8 @@ class CubicHorizontalKernel {
     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);
+
+        make_valid_stencil(x, y, ws.stencil);
         const auto& wj = ws.weights.weights_j;
 
         idx_t pos = row * stencil_size();