Merge branch 'develop' into feature/ecmwf-ukv-yaml

* develop: (21 commits)
  Add new yml file for variable resolution projection To run the update_uid.py under linux I should add the option shell=True in Popen
  Make VortexRollup function publicly available, and remove it from various places (#87)
  More testing and implementations for VariableResolutionProjection (#89)
  Fixup PointXYZ operator /= to be bitreproducible with PointXYZ::div
  Add *= and /= operators to PointXYZ with scalar
  ATLAS-345 Use new eckit::linalg::LinearAlgebraSparse and eckit::linalg::LinearAlgebraDense API (see ECKIT-585)
  Prevent optimisation with NVHPC compiler in VariableResolutionProjection
  Adapt atlas_test_field_missingvalue as there is no standard on comparisons with NaN
  Workaround for Cray 8.7 compiler bug resulting in FE_DIVBYZERO in atlas_test_jacobian in Release build
  ATLAS-327 Workaround for Cray 8.7 compiler bug with CXX flags "-O2 -hfp1 -g -DNDEBUG"
  Fixup previous commit
  Allow use of build_cells_remote_index() when halo was already built
  Apply clang-format
  Addressed reviewer comments.
  Permitted variable sized halos on functionspaces.
  Added ghost field to CubedSphereCellColumns functionspace.
  ATLAS-344 Fix atlas_test_healpixmeshgen CI with nvhpc-21.9
  ATLAS-344 Fix MatchingMeshPartitionerLonLatPolygon to work with non-matching domains
  ATLAS-344 PolygonXY with Point2 member variables
  ATLAS-344 Printing of util::PolygonCoordinates
  ...

src/apps/atlas-meshgen.cc

@@ -125,6 +125,7 @@ class Meshgen2Gmsh : public AtlasTool {
     std::string key;
     long halo;
     bool edges;
+    bool cells;
     bool brick;
     bool stats;
     bool info;
@@ -165,6 +166,7 @@ Meshgen2Gmsh::Meshgen2Gmsh(int argc, char** argv): AtlasTool(argc, argv) {
     add_option(new Separator("Advanced"));
     add_option(new SimpleOption<long>("halo", "Halo size"));
     add_option(new SimpleOption<bool>("edges", "Build edge datastructure"));
+    add_option(new SimpleOption<bool>("cells", "Build cells datastructure"));
     add_option(new SimpleOption<bool>("brick", "Build brick dual mesh"));
     add_option(new SimpleOption<bool>("stats", "Write statistics file"));
     add_option(new SimpleOption<bool>("info", "Write Info"));
@@ -203,6 +205,8 @@ int Meshgen2Gmsh::execute(const Args& args) {
 
     edges = false;
     args.get("edges", edges);
+    cells = false;
+    args.get("cells", cells);
     stats = false;
     args.get("stats", stats);
     info = false;
@@ -280,6 +284,7 @@ int Meshgen2Gmsh::execute(const Args& args) {
         throw;
     }
 
+
     if (grid.projection().units() == "degrees") {
         functionspace::NodeColumns nodes_fs(mesh, option::halo(halo));
     }
@@ -299,6 +304,10 @@ int Meshgen2Gmsh::execute(const Args& args) {
         }
     }
 
+    if (cells) {
+        functionspace::CellColumns cells_fs(mesh, option::halo(halo));
+    }
+
     if (stats) {
         build_statistics(mesh);
     }

src/atlas/CMakeLists.txt

@@ -728,6 +728,8 @@ util/detail/BlackMagic.h
 util/detail/Cache.h
 util/detail/Debug.h
 util/detail/KDTree.h
+util/function/VortexRollup.h
+util/function/VortexRollup.cc
 )
 
 list( APPEND atlas_internals_srcs

src/atlas/functionspace/CellColumns.cc

@@ -272,10 +272,10 @@ CellColumns::CellColumns(const Mesh& mesh, const eckit::Configuration& config):
     mesh::actions::build_cells_parallel_fields(mesh_);
     mesh::actions::build_periodic_boundaries(mesh_);
 
-    if (halo_.size() > 0) {
-        mesh::actions::build_halo(mesh_, halo_.size());
-        nb_cells_ = get_nb_cells_from_metadata();
-    }
+
+    mesh::actions::build_halo(mesh_, halo_.size());
+    nb_cells_ = get_nb_cells_from_metadata();
+
     if (!nb_cells_) {
         nb_cells_ = mesh.cells().size();
     }
@@ -575,6 +575,10 @@ Field CellColumns::global_index() const {
     return mesh_.cells().global_index();
 }
 
+Field CellColumns::ghost() const {
+    return mesh_.cells().field("ghost");
+}
+
 //------------------------------------------------------------------------------
 //------------------------------------------------------------------------------
 //------------------------------------------------------------------------------

src/atlas/functionspace/CellColumns.h

@@ -93,6 +93,8 @@ class CellColumns : public functionspace::FunctionSpaceImpl {
 
     Field global_index() const override;
 
+    Field ghost() const override;
+
 private:  // methods
     idx_t config_size(const eckit::Configuration& config) const;
     array::DataType config_datatype(const eckit::Configuration&) const;

src/atlas/functionspace/CubedSphereColumns.cc

@@ -65,7 +65,7 @@ class CubedSphereStructureCache : public util::Cache<std::string, detail::CubedS
         ATLAS_ASSERT(mesh);
         auto& mesh_impl = *mesh.get();
         registerMesh(mesh_impl);
-        creator_type creator                 = std::bind(&CubedSphereStructureCache::create, mesh);
+        creator_type creator = std::bind(&CubedSphereStructureCache::create, mesh, functionspace->size());
         util::ObjectHandle<value_type> value = Base::get_or_create(key(mesh_impl), creator);
         return value;
     }
@@ -78,8 +78,8 @@ class CubedSphereStructureCache : public util::Cache<std::string, detail::CubedS
         return key.str();
     }
 
-    static value_type* create(const Mesh& mesh) {
-        value_type* value = new value_type(getTij<BaseFunctionSpace>(mesh), getGhost<BaseFunctionSpace>(mesh));
+    static value_type* create(const Mesh& mesh, idx_t size) {
+        value_type* value = new value_type(getTij<BaseFunctionSpace>(mesh), getGhost<BaseFunctionSpace>(mesh), size);
         return value;
     }
 };
@@ -121,13 +121,8 @@ idx_t CubedSphereColumns<BaseFunctionSpace>::invalid_index() const {
 }
 
 template <typename BaseFunctionSpace>
-idx_t CubedSphereColumns<BaseFunctionSpace>::nb_elems() const {
-    return cubedSphereColumnsHandle_.get()->nb_elems();
-}
-
-template <typename BaseFunctionSpace>
-idx_t CubedSphereColumns<BaseFunctionSpace>::nb_owned_elems() const {
-    return cubedSphereColumnsHandle_.get()->nb_owned_elems();
+idx_t CubedSphereColumns<BaseFunctionSpace>::sizeOwned() const {
+    return cubedSphereColumnsHandle_.get()->sizeOwned();
 }
 
 template <typename BaseFunctionSpace>
@@ -160,11 +155,6 @@ Field CubedSphereColumns<BaseFunctionSpace>::tij() const {
     return cubedSphereColumnsHandle_.get()->tij();
 }
 
-template <typename BaseFunctionSpace>
-Field CubedSphereColumns<BaseFunctionSpace>::ghost() const {
-    return cubedSphereColumnsHandle_.get()->ghost();
-}
-
 // Explicit instantiation of template classes.
 template class CubedSphereColumns<CellColumns>;
 template class CubedSphereColumns<NodeColumns>;

src/atlas/functionspace/CubedSphereColumns.h

@@ -38,11 +38,8 @@ class CubedSphereColumns : public BaseFunctionSpace {
     /// Invalid index.
     idx_t invalid_index() const;
 
-    /// Get total number of elements.
-    idx_t nb_elems() const;
-
     /// Get number of owned elements.
-    idx_t nb_owned_elems() const;
+    idx_t sizeOwned() const;
 
     /// i lower bound for tile t (including halo)
     idx_t i_begin(idx_t t) const;
@@ -60,16 +57,12 @@ class CubedSphereColumns : public BaseFunctionSpace {
     /// Return tij field.
     Field tij() const;
 
-    /// Return ghost field.
-    Field ghost() const;
-
 private:
     class For {
     public:
         For(const CubedSphereColumns<BaseFunctionSpace>& functionSpace, const util::Config& config = util::NoConfig()):
             functionSpace_{functionSpace},
-            indexMax_{config.getBool("include_halo", false) ? functionSpace.nb_elems()
-                                                            : functionSpace.nb_owned_elems()},
+            indexMax_{config.getBool("include_halo", false) ? functionSpace.size() : functionSpace.sizeOwned()},
             levels_{config.getInt("levels", functionSpace_.levels())},
             tijView_(array::make_view<idx_t, 2>(functionSpace_.tij())) {}
 
@@ -78,7 +71,7 @@ class CubedSphereColumns : public BaseFunctionSpace {
         using EnableFunctor =
             typename std::enable_if<std::is_convertible<FuncType, std::function<void(ArgTypes...)>>::value>::type*;
 
-        // Functor: void f( index, t, i, j, k)
+        // Functor: void f(index, t, i, j, k)
         template <typename Functor, EnableFunctor<Functor, idx_t, idx_t, idx_t, idx_t, idx_t> = nullptr>
         void operator()(const Functor& f) const {
             using namespace meshgenerator::detail::cubedsphere;
@@ -94,7 +87,7 @@ class CubedSphereColumns : public BaseFunctionSpace {
             }
         }
 
-        // Functor: void f( index, t, i, j)
+        // Functor: void f(index, t, i, j)
         template <typename Functor, EnableFunctor<Functor, idx_t, idx_t, idx_t, idx_t> = nullptr>
         void operator()(const Functor& f) const {
             using namespace meshgenerator::detail::cubedsphere;
@@ -108,7 +101,7 @@ class CubedSphereColumns : public BaseFunctionSpace {
             }
         }
 
-        // Functor: void f( index, k)
+        // Functor: void f(index, k)
         template <typename Functor, EnableFunctor<Functor, idx_t, idx_t> = nullptr>
         void operator()(const Functor& f) const {
             using namespace meshgenerator::detail::cubedsphere;
@@ -121,7 +114,7 @@ class CubedSphereColumns : public BaseFunctionSpace {
             }
         }
 
-        // Functor: void f( index )
+        // Functor: void f(index )
         template <typename Functor, EnableFunctor<Functor, idx_t> = nullptr>
         void operator()(const Functor& f) const {
             using namespace meshgenerator::detail::cubedsphere;

src/atlas/functionspace/detail/CubedSphereStructure.cc

@@ -32,16 +32,15 @@ CubedSphereStructure::BoundingBox::BoundingBox() {
     jEnd   = std::numeric_limits<idx_t>::min();
 }
 
-CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost): tij_(tij), ghost_(ghost) {
+CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost, idx_t size):
+    tij_(tij), ghost_(ghost), nElems_(size) {
     ATLAS_TRACE();
     Log::debug() << "CubedSphereStructure bounds checking is set to " + std::to_string(checkBounds) << std::endl;
 
     // Make array views.
     const auto tijView_   = array::make_view<idx_t, 2>(tij_);
     const auto ghostView_ = array::make_view<int, 1>(ghost_);
 
-    nElems_ = tijView_.shape(0);
-
     // loop over tij and find min and max ij bounds.
     for (idx_t index = 0; index < nElems_; ++index) {
         const size_t t = static_cast<size_t>(tijView_(index, Coordinates::T));
@@ -76,11 +75,11 @@ CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost)
     }
 }
 
-idx_t CubedSphereStructure::nb_elems() const {
+idx_t CubedSphereStructure::size() const {
     return nElems_;
 }
 
-idx_t CubedSphereStructure::nb_owned_elems() const {
+idx_t CubedSphereStructure::sizeOwned() const {
     return nOwnedElems_;
 }
 

src/atlas/functionspace/detail/CubedSphereStructure.h

@@ -27,16 +27,16 @@ namespace detail {
 class CubedSphereStructure : public util::Object {
 public:
     CubedSphereStructure() = default;
-    CubedSphereStructure(const Field& tij, const Field& ghost);
+    CubedSphereStructure(const Field& tij, const Field& ghost, idx_t size);
 
     /// Invalid index.
     static constexpr idx_t invalid_index() { return -1; }
 
     /// Number of elements.
-    idx_t nb_elems() const;
+    idx_t size() const;
 
     /// Number of owned elements.
-    idx_t nb_owned_elems() const;
+    idx_t sizeOwned() const;
 
     /// i lower bound for tile t (including halo)
     idx_t i_begin(idx_t) const;

src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.cc

@@ -269,6 +269,10 @@ void cap_colats(int N, int n_collars, const double& c_polar, int n_regions[], do
     c_caps[n_collars + 1] = M_PI;
 }
 
+// Disable optimisation for Cray (See ATLAS-327)
+#ifdef _CRAYC
+#pragma _CRI noopt
+#endif
 void eq_caps(int N, std::vector<int>& n_regions, std::vector<double>& s_cap) {
     //
     // eq_regions uses the zonal equal area sphere partitioning algorithm to
@@ -291,14 +295,14 @@ void eq_caps(int N, std::vector<int>& n_regions, std::vector<double>& s_cap) {
         // Given N, determine c_polar
         // the colatitude of the North polar spherical cap.
         //
-        double c_polar = polar_colat(N);
+        const double c_polar = polar_colat(N);
 
         //
         // Given N, determine the ideal angle for spherical collars.
         // Based on N, this ideal angle, and c_polar,
         // determine n_collars, the number of collars between the polar caps.
         //
-        int n_collars = num_collars(N, c_polar, ideal_collar_angle(N));
+        const int n_collars = num_collars(N, c_polar, ideal_collar_angle(N));
 
         // int n_regions_ns=n_collars+2;
         //
@@ -337,6 +341,10 @@ void eq_caps(int N, std::vector<int>& n_regions, std::vector<double>& s_cap) {
     }
     // int n_regions_ew=maxval(n_regions(:));
 }
+// Reenable optimisation for Cray (See ATLAS-327)
+#ifdef _CRAYC
+#pragma _CRI opt
+#endif
 
 void eq_regions(int N, double xmin[], double xmax[], double ymin[], double ymax[]) {
     // EQ_REGIONS Recursive zonal equal area (EQ) partition of sphere

src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc

@@ -10,6 +10,8 @@
 
 #include "atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h"
 
+#include <iomanip>
+#include <sstream>
 #include <vector>
 
 #include "eckit/config/Resource.h"
@@ -19,6 +21,7 @@
 #include "atlas/grid/Iterator.h"
 #include "atlas/mesh/Nodes.h"
 #include "atlas/parallel/mpi/mpi.h"
+#include "atlas/parallel/omp/fill.h"
 #include "atlas/runtime/Exception.h"
 #include "atlas/runtime/Log.h"
 #include "atlas/util/CoordinateEnums.h"
@@ -44,36 +47,63 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti
 
     Log::debug() << "MatchingMeshPartitionerLonLatPolygon::partition" << std::endl;
 
-    // FIXME: THIS IS A HACK! the coordinates include North/South Pole (first/last
-    // partitions only)
-    bool includesNorthPole = (mpi_rank == 0);
-    bool includesSouthPole = (mpi_rank == mpi_size - 1);
-
     const util::PolygonXY poly{prePartitionedMesh_.polygon(0)};
+
+    double west = poly.coordinatesMin().x();
+    double east = poly.coordinatesMax().x();
+    comm.allReduceInPlace(west, eckit::mpi::Operation::MIN);
+    comm.allReduceInPlace(east, eckit::mpi::Operation::MAX);
+
     Projection projection = prePartitionedMesh_.projection();
+    omp::fill(partitioning, partitioning + grid.size(), -1);
 
-    {
-        eckit::ProgressTimer timer("Partitioning", grid.size(), "point", double(10), atlas::Log::trace());
+    auto compute = [&](double west) {
         size_t i = 0;
 
         for (PointLonLat P : grid.lonlat()) {
-            ++timer;
-            projection.lonlat2xy(P);
-            const bool atThePole = (includesNorthPole && P[LAT] >= poly.coordinatesMax()[LAT]) ||
-                                   (includesSouthPole && P[LAT] < poly.coordinatesMin()[LAT]);
-
-            partitioning[i++] = atThePole || poly.contains(P) ? mpi_rank : -1;
+            if (partitioning[i] < 0) {
+                projection.lonlat2xy(P);
+                P.normalise(west);
+                partitioning[i] = poly.contains(P) ? mpi_rank : -1;
+            }
+            ++i;
         }
-    }
+        // Synchronize partitioning
+        comm.allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX);
 
-    // Synchronize partitioning, do a sanity check
-    comm.allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX);
-    const int min = *std::min_element(partitioning, partitioning + grid.size());
+        return *std::min_element(partitioning, partitioning + grid.size());
+    };
+
+    int min              = compute(east - 360.);
+    constexpr double eps = 1.e-10;
+    if (min < 0 && east - west > 360. + eps) {
+        min = compute(west - eps);
+    }
     if (min < 0) {
-        throw_Exception(
-            "Could not find partition for target node (source "
-            "mesh does not contain all target grid points)",
-            Here());
+        size_t i            = 0;
+        size_t max_failures = grid.size();
+        std::vector<size_t> failed_index;
+        std::vector<PointLonLat> failed_lonlat;
+        failed_index.reserve(max_failures);
+        failed_lonlat.reserve(max_failures);
+        for (PointLonLat P : grid.lonlat()) {
+            if (partitioning[i] < 0) {
+                failed_index.emplace_back(i);
+                failed_lonlat.emplace_back(P);
+            }
+            ++i;
+        }
+        size_t nb_failures = failed_index.size();
+        std::stringstream err;
+        err << "Could not find partition of " << nb_failures
+            << " target grid points (source "
+               "mesh does not contain all target grid points)\n"
+               "Failed target grid points with global index:\n";
+        for (size_t n = 0; n < nb_failures; ++n) {
+            err << "  - " << std::setw(10) << std::left << failed_index[n] + 1 << " {lon,lat} : " << failed_lonlat[n]
+                << "\n";
+        }
+        throw_Exception(err.str(), Here());
     }
 }