Merge branch 'develop' into feature/parallel_transport

CHANGELOG.md

@@ -7,7 +7,21 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 
 ## [Unreleased]
 
-## [0.35.1] - 2023-24-10
+## [0.36.0] - 2023-12-11
+### Added
+- Add TriangularMeshBuilder with Fortran API, so far for serial meshes only
+- Add Fortran API for CellCollumns functionspace
+- Add EqualAreaPartitioner which is geometry based rather than loadbalanced like EqualRegionsPartitioner
+
+### Fixed
+- Compilation with Intel LLVM compiler
+- Fix 180 degrees phase shift error in MDPI_gulfstream function
+
+### Changed
+- Update install scripts
+- Preparation for using eckit::codec as backend for atlas_io
+
+## [0.35.1] - 2023-10-24
 ### Added
 - Don't output with output::Gmsh the triangle elements with wrong orientation when coordinates are "lonlat"
 - Add control to skip Gmsh-output of triangles with too large edge length ratio
@@ -500,6 +514,7 @@ Fix StructuredInterpolation2D with retry for failed stencils
 ## 0.13.0 - 2018-02-16
 
 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
+[0.36.0]: https://github.com/ecmwf/atlas/compare/0.35.1...0.36.0
 [0.35.1]: https://github.com/ecmwf/atlas/compare/0.35.0...0.35.1
 [0.35.0]: https://github.com/ecmwf/atlas/compare/0.34.0...0.35.0
 [0.34.0]: https://github.com/ecmwf/atlas/compare/0.33.0...0.34.0

VERSION

@@ -1 +1 @@
-0.35.1
+0.36.0

src/apps/atlas-grid-points.cc

@@ -55,7 +55,7 @@ Program::Program(int argc, char** argv): AtlasTool(argc, argv) {
     add_option(new SimpleOption<std::string>("field","Field to output. [\"lonlat\"(D),\"index\",\"partition\"]"));
     add_option(new Separator("Advanced"));
     add_option(new SimpleOption<std::string>("partitioner.type",
-        "Partitioner [equal_regions,checkerboard,equal_bands,regular_bands]"));
+        "Partitioner [equal_regions,equal_area,checkerboard,equal_bands,regular_bands]"));
     add_option(new SimpleOption<long>("partition", "partition [0:partitions-1]"));
     add_option(new SimpleOption<long>("partitions", "Number of partitions"));
     add_option(new SimpleOption<long>("json.precision", "Number of digits after decimal in output"));

src/atlas/CMakeLists.txt

@@ -288,6 +288,8 @@ grid/detail/partitioner/CubedSpherePartitioner.cc
 grid/detail/partitioner/CubedSpherePartitioner.h
 grid/detail/partitioner/EqualBandsPartitioner.cc
 grid/detail/partitioner/EqualBandsPartitioner.h
+grid/detail/partitioner/EqualAreaPartitioner.cc
+grid/detail/partitioner/EqualAreaPartitioner.h
 grid/detail/partitioner/EqualRegionsPartitioner.cc
 grid/detail/partitioner/EqualRegionsPartitioner.h
 grid/detail/partitioner/MatchingMeshPartitioner.h

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

@@ -0,0 +1,77 @@
+/*
+ * (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/grid/detail/partitioner/EqualAreaPartitioner.h"
+
+#include "atlas/grid.h"
+#include "atlas/grid/Iterator.h"
+#include "atlas/util/Constants.h"
+
+namespace atlas {
+namespace grid {
+namespace detail {
+namespace partitioner {
+
+EqualAreaPartitioner::EqualAreaPartitioner():
+    Partitioner(), partitioner_() {
+}
+
+EqualAreaPartitioner::EqualAreaPartitioner(int N):
+    EqualAreaPartitioner(N, util::NoConfig()) {
+}
+
+EqualAreaPartitioner::EqualAreaPartitioner(int N, const eckit::Parametrisation& config):
+    Partitioner(N,config), partitioner_(N,config) {
+}
+
+EqualAreaPartitioner::EqualAreaPartitioner(const eckit::Parametrisation& config):
+    Partitioner(config), partitioner_(config) {
+}
+
+void EqualAreaPartitioner::partition(const Grid& grid, int part[]) const {
+
+    if( partitioner_.coordinates_ == EqualRegionsPartitioner::Coordinates::XY && StructuredGrid(grid) ) {
+        StructuredGrid g(grid);
+        size_t n = 0;
+        for (idx_t j=0; j<g.ny(); ++j) {
+            const double lat = g.y(j) * util::Constants::degreesToRadians();
+            int b = partitioner_.band(lat);
+            int p = 0;
+            for (int k = 0; k < b; ++k) {
+                p += partitioner_.nb_regions(k);
+            }
+            idx_t nx = g.nx(j);
+            for (idx_t i=0; i<nx; ++i) {
+                const double lon = g.x(i,j) * util::Constants::degreesToRadians();
+                part[n++] = p + partitioner_.sector(b, lon);
+            }
+        }
+    }
+    else {
+        size_t j{0};
+        for (PointLonLat p : grid.lonlat()) {
+            p.lon() *= util::Constants::degreesToRadians();
+            p.lat() *= util::Constants::degreesToRadians();
+            part[j++] = partitioner_.partition(p.lon(), p.lat());
+        }
+    }
+
+}
+
+
+}  // namespace partitioner
+}  // namespace detail
+}  // namespace grid
+}  // namespace atlas
+
+namespace {
+atlas::grid::detail::partitioner::PartitionerBuilder<atlas::grid::detail::partitioner::EqualAreaPartitioner>
+    __EqualArea("equal_area");
+}

src/atlas/grid/detail/partitioner/EqualAreaPartitioner.h

@@ -0,0 +1,43 @@
+/*
+ * (C) Copyright 2023 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
+
+#include <vector>
+
+#include "atlas/grid/detail/partitioner/EqualRegionsPartitioner.h"
+#include "atlas/grid/detail/partitioner/Partitioner.h"
+
+namespace atlas {
+namespace grid {
+namespace detail {
+namespace partitioner {
+
+class EqualAreaPartitioner : public Partitioner {
+public:
+    EqualAreaPartitioner();
+    EqualAreaPartitioner(int N);
+    EqualAreaPartitioner(int N, const eckit::Parametrisation& config);
+    EqualAreaPartitioner(const eckit::Parametrisation& config);
+
+    using Partitioner::partition;
+    void partition(const Grid&, int part[]) const override;
+
+    std::string type() const override { return "equal_area"; }
+
+private:
+    EqualRegionsPartitioner partitioner_;
+};
+
+
+}  // namespace partitioner
+}  // namespace detail
+}  // namespace grid
+}  // namespace atlas

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

@@ -495,14 +495,16 @@ int EqualRegionsPartitioner::band(const double& y) const {
 }
 
 int EqualRegionsPartitioner::sector(int band, const double& x) const {
+    constexpr double two_pi = 2. * M_PI;
+    constexpr double div_two_pi = 1. / (2. * M_PI + 1.e-8);
     double xreg = x;
     if (x < 0.) {
-        xreg += 2. * M_PI;
+        xreg += two_pi;
     }
-    else if (x > 2. * M_PI) {
-        xreg -= 2. * M_PI;
+    else if (x > two_pi) {
+        xreg -= two_pi;
     }
-    return std::floor(xreg * sectors_[band] / (2. * M_PI + 1e-8));
+    return std::floor(xreg * sectors_[band] * div_two_pi);
 }
 
 void EqualRegionsPartitioner::where(int partition, int& band, int& sector) const {

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

@@ -75,6 +75,7 @@ namespace partitioner {
 void eq_caps(int N, std::vector<int>& n_regions, std::vector<double>& s_cap);
 void eq_regions(int N, double xmin[], double xmax[], double ymin[], double ymax[]);
 
+class EqualAreaPartitioner;
 class EqualRegionsPartitioner : public Partitioner {
 public:
     EqualRegionsPartitioner();
@@ -121,15 +122,16 @@ class EqualRegionsPartitioner : public Partitioner {
     // algorithm is used internally
     void partition(int nb_nodes, NodeInt nodes[], int part[]) const;
 
-    // x and y in radians
-    int partition(const double& x, const double& y) const;
-
+    friend class EqualAreaPartitioner;
     // y in radians
     int band(const double& y) const;
 
     // x in radians
     int sector(int band, const double& x) const;
 
+    // x and y in radians
+    int partition(const double& x, const double& y) const;
+
 private:
     int N_;
     std::vector<double> bands_;

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

@@ -16,7 +16,7 @@
 
 #include "atlas/array/ArrayView.h"
 #include "atlas/field/Field.h"
-#include "atlas/grid/Grid.h"
+#include "atlas/grid.h"
 #include "atlas/mesh/Elements.h"
 #include "atlas/mesh/Mesh.h"
 #include "atlas/mesh/Nodes.h"
@@ -93,21 +93,23 @@ void MatchingMeshPartitionerBruteForce::partition(const Grid& grid, int partitio
     auto lonlat_src = array::make_view<double, 2>(prePartitionedMesh_.nodes().lonlat());
 
     std::vector<PointLonLat> coordinates;
+    coordinates.reserve(lonlat_src.shape(0));
     PointLonLat coordinatesMin = PointLonLat(lonlat_src(0, LON), lonlat_src(0, LAT));
     PointLonLat coordinatesMax = coordinatesMin;
 
-    for (idx_t i = 0; i < lonlat_src.size(); ++i) {
+    for (idx_t i = 0; i < lonlat_src.shape(0); ++i) {
         PointLonLat A(lonlat_src(i, LON), lonlat_src(i, LAT));
         coordinatesMin = PointLonLat::componentsMin(coordinatesMin, A);
         coordinatesMax = PointLonLat::componentsMax(coordinatesMax, A);
-        coordinates.push_back(A);
+        coordinates.emplace_back(A);
     }
 
     {
         eckit::ProgressTimer timer("Partitioning target", grid.size(), "point", double(10), atlas::Log::trace());
-        for (idx_t i = 0; i < grid.size(); ++i, ++timer) {
+        auto grid_iter = grid.lonlat().begin();
+        for (idx_t i = 0; i < grid.size(); ++i, ++grid_iter) {
             partitioning[i] = -1;
-            const PointLonLat& P(coordinates[i]);
+            const PointLonLat& P = *grid_iter;
 
             if (coordinatesMin[LON] <= P[LON] && P[LON] <= coordinatesMax[LON] && coordinatesMin[LAT] <= P[LAT] &&
                 P[LAT] <= coordinatesMax[LAT]) {

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

@@ -28,6 +28,9 @@
 #include "atlas/util/CoordinateEnums.h"
 #include "atlas/util/PolygonXY.h"
 
+
+#include "atlas/util/KDTree.h"
+
 namespace atlas {
 namespace grid {
 namespace detail {
@@ -37,6 +40,12 @@ namespace {
 PartitionerBuilder<MatchingMeshPartitionerLonLatPolygon> __builder("lonlat-polygon");
 }
 
+MatchingMeshPartitionerLonLatPolygon::MatchingMeshPartitionerLonLatPolygon(const Mesh& mesh, const eckit::Parametrisation& config):
+    MatchingMeshPartitioner(mesh, config) {
+        config.get("fallback_nearest", fallback_nearest_);
+    }
+
+
 void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int partitioning[]) const {
     const auto& comm   = mpi::comm(prePartitionedMesh_.mpi_comm());
     const int mpi_rank = int(comm.rank());
@@ -122,6 +131,7 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti
         }
         return false;
     }();
+
     if (min < 0) {
         size_t i            = 0;
         size_t max_failures = grid.size();
@@ -137,20 +147,41 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti
             ++i;
         }
         size_t nb_failures = failed_index.size();
-        std::stringstream err;
-        err.precision(20);
-        err << "Could not find partition of " << nb_failures
-            << " target grid points (source mesh does not contain all target grid points)\n"
-            << "Tried first normalizing coordinates with west=" << east - 360.;
-        if (second_try) {
-            err << "Tried second time normalizing coordinates with west=" << west - eps << "\n";
+
+        if (fallback_nearest_) {
+            util::IndexKDTree3D kdtree;
+            kdtree.reserve(grid.size());
+            size_t j=0;
+            for (auto& p: grid.lonlat()) {
+                if (partitioning[j] >= 0) {
+                    kdtree.insert(p,partitioning[j]);
+                }
+                ++j;
+            }
+            kdtree.build();
+            for (size_t n = 0; n < nb_failures; ++n) {
+                auto closest = kdtree.closestPoint(failed_lonlat[n]);
+                partitioning[failed_index[n]] = closest.payload();
+            }
         }
-        err << "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";
+        else {
+            std::stringstream err;
+            err.precision(20);
+            err << "Could not find partition of " << nb_failures
+                << " target grid points (source mesh does not contain all target grid points)\n"
+                << "Tried first normalizing coordinates with west=" << east - 360. << "\n";
+            if (second_try) {
+                err << "Tried second time normalizing coordinates with west=" << west - eps << "\n";
+            }
+
+            err << "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";
+            }
+            // prePartitionedMesh_.polygon(0).outputPythonScript("partitions.py");
+            throw_Exception(err.str(), Here());
         }
-        throw_Exception(err.str(), Here());
     }
 }
 

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

@@ -22,7 +22,7 @@ class MatchingMeshPartitionerLonLatPolygon : public MatchingMeshPartitioner {
     static std::string static_type() { return "lonlat-polygon"; }
 
 public:
-    MatchingMeshPartitionerLonLatPolygon(const Mesh& mesh, const eckit::Parametrisation& config): MatchingMeshPartitioner(mesh, config) {}
+    MatchingMeshPartitionerLonLatPolygon(const Mesh& mesh, const eckit::Parametrisation& config);
 
 
     MatchingMeshPartitionerLonLatPolygon(): MatchingMeshPartitioner() {}
@@ -44,6 +44,9 @@ class MatchingMeshPartitionerLonLatPolygon : public MatchingMeshPartitioner {
     void partition(const Grid& grid, int partitioning[]) const;
 
     virtual std::string type() const { return static_type(); }
+
+private:
+    bool fallback_nearest_{false};
 };
 
 }  // namespace partitioner

src/atlas/util/GridPointsJSONWriter.cc

@@ -98,6 +98,11 @@ void GridPointsJSONWriter::write(std::ostream& out, eckit::Channel& info) const
 //------------------------------------------------------------------------------
 
 void GridPointsJSONWriter::write(std::ostream& out, std::ostream* info) const {
+
+    if (field_ == "none") {
+        return;
+    }
+
     int points_newline = pretty_ ? true : false;
     int points_indent = pretty_ ? 2 : 0;
     int partition_indent = 0;

src/atlas/util/function/MDPI_functions.cc

@@ -93,7 +93,7 @@ double MDPI_gulfstream(double lon, double 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.;
+        double gf_per_lon = lon;
         while (gf_per_lon > 180.) {
             gf_per_lon -= 360.;
         }

src/atlas_f/atlas_module.F90

@@ -151,7 +151,8 @@ module atlas_module
 use atlas_trace_module, only : &
     & atlas_Trace
 use atlas_functions_module
-
+use atlas_Redistribution_module, only : &
+    & atlas_Redistribution
 use fckit_log_module,  only: atlas_log => fckit_log
 
 implicit none