From f3f83fcefbf7a6ca7c7fb62b28f67de66ebec3f0 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 6 Dec 2023 18:08:15 +0000 Subject: [PATCH 1/8] Fix MDPI_gulfstream: 180 degrees phase shift corrected --- src/atlas/util/function/MDPI_functions.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/util/function/MDPI_functions.cc b/src/atlas/util/function/MDPI_functions.cc index f570ac5b7..58713a037 100644 --- a/src/atlas/util/function/MDPI_functions.cc +++ b/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.; } From 68b09acef60294f10e455dfe905a86a7c96c9620 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 6 Dec 2023 17:02:20 +0000 Subject: [PATCH 2/8] Fix MatchingMeshPartitionerBruteForce --- .../partitioner/MatchingMeshPartitionerBruteForce.cc | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerBruteForce.cc b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerBruteForce.cc index c99afe0c1..408014ba8 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerBruteForce.cc +++ b/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(prePartitionedMesh_.nodes().lonlat()); std::vector 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]) { From de9c9fde6d17a6e18c17a4c5176f2da10ad33f9f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 7 Dec 2023 11:05:17 +0000 Subject: [PATCH 3/8] Add fallback mechanism to MatchingMeshPartitionerLonLatPolygon --- .../MatchingMeshPartitionerLonLatPolygon.cc | 55 +++++++++++++++---- .../MatchingMeshPartitionerLonLatPolygon.h | 5 +- 2 files changed, 47 insertions(+), 13 deletions(-) diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc index d517c8b2b..0e06d587e 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc +++ b/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 __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()); } } diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h index f0b11e5d3..4dd45f0ec 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h +++ b/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 From bdc4ccf276f8bb944251ee03fa95d70323f08117 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 7 Dec 2023 12:02:08 +0000 Subject: [PATCH 4/8] Add EqualAreaPartitioner which is geometry based compared to EqualRegionsPartitioner which is loadbalanced --- src/atlas/CMakeLists.txt | 2 + .../partitioner/EqualAreaPartitioner.cc | 56 +++++++++++++++++++ .../detail/partitioner/EqualAreaPartitioner.h | 43 ++++++++++++++ .../partitioner/EqualRegionsPartitioner.h | 8 ++- 4 files changed, 106 insertions(+), 3 deletions(-) create mode 100644 src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc create mode 100644 src/atlas/grid/detail/partitioner/EqualAreaPartitioner.h diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 3627dd60d..02d2ec38d 100644 --- a/src/atlas/CMakeLists.txt +++ b/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 diff --git a/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc b/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc new file mode 100644 index 000000000..ced5a3957 --- /dev/null +++ b/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc @@ -0,0 +1,56 @@ +/* + * (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/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& g, int part[]) const { + size_t j{0}; + for (PointLonLat p : g.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 + __EqualArea("equal_area"); +} diff --git a/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.h b/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.h new file mode 100644 index 000000000..87f43e024 --- /dev/null +++ b/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 + +#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 diff --git a/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.h b/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.h index a795a6f15..511879993 100644 --- a/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.h +++ b/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.h @@ -75,6 +75,7 @@ namespace partitioner { void eq_caps(int N, std::vector& n_regions, std::vector& 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 bands_; From 5517a960965886161925d7a8a2482433605107c5 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 7 Dec 2023 15:15:51 +0000 Subject: [PATCH 5/8] Performance improvement in EqualAreaPartitioner for StructuredGrid --- src/apps/atlas-grid-points.cc | 2 +- .../partitioner/EqualAreaPartitioner.cc | 35 +++++++++++++++---- .../partitioner/EqualRegionsPartitioner.cc | 10 +++--- src/atlas/util/GridPointsJSONWriter.cc | 5 +++ 4 files changed, 40 insertions(+), 12 deletions(-) diff --git a/src/apps/atlas-grid-points.cc b/src/apps/atlas-grid-points.cc index 86cf4c674..8a74c83c0 100644 --- a/src/apps/atlas-grid-points.cc +++ b/src/apps/atlas-grid-points.cc @@ -55,7 +55,7 @@ Program::Program(int argc, char** argv): AtlasTool(argc, argv) { add_option(new SimpleOption("field","Field to output. [\"lonlat\"(D),\"index\",\"partition\"]")); add_option(new Separator("Advanced")); add_option(new SimpleOption("partitioner.type", - "Partitioner [equal_regions,checkerboard,equal_bands,regular_bands]")); + "Partitioner [equal_regions,equal_area,checkerboard,equal_bands,regular_bands]")); add_option(new SimpleOption("partition", "partition [0:partitions-1]")); add_option(new SimpleOption("partitions", "Number of partitions")); add_option(new SimpleOption("json.precision", "Number of digits after decimal in output")); diff --git a/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc b/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc index ced5a3957..ac585e316 100644 --- a/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc +++ b/src/atlas/grid/detail/partitioner/EqualAreaPartitioner.cc @@ -10,7 +10,7 @@ #include "atlas/grid/detail/partitioner/EqualAreaPartitioner.h" -#include "atlas/grid/Grid.h" +#include "atlas/grid.h" #include "atlas/grid/Iterator.h" #include "atlas/util/Constants.h" @@ -35,13 +35,34 @@ EqualAreaPartitioner::EqualAreaPartitioner(const eckit::Parametrisation& config) Partitioner(config), partitioner_(config) { } -void EqualAreaPartitioner::partition(const Grid& g, int part[]) const { - size_t j{0}; - for (PointLonLat p : g.lonlat()) { - p.lon() *= util::Constants::degreesToRadians(); - p.lat() *= util::Constants::degreesToRadians(); - part[j++] = partitioner_.partition(p.lon(), p.lat()); +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 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 { diff --git a/src/atlas/util/GridPointsJSONWriter.cc b/src/atlas/util/GridPointsJSONWriter.cc index e7066ad77..7888a0c9a 100644 --- a/src/atlas/util/GridPointsJSONWriter.cc +++ b/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; From 2ffed1dd47a5361a7dc058da6de6c62387b4fa86 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 7 Dec 2023 15:32:18 +0000 Subject: [PATCH 6/8] Add atlas_Redistribution to atlas_module --- src/atlas_f/atlas_module.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/atlas_f/atlas_module.F90 b/src/atlas_f/atlas_module.F90 index 975ea4e92..97ddb45ae 100644 --- a/src/atlas_f/atlas_module.F90 +++ b/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 From f2e0bc286f23e13736796f58d67abe98f882771d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 6 Dec 2023 15:18:50 +0000 Subject: [PATCH 7/8] Version 0.36.0 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 731b95d7f..93d4c1ef0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.35.1 +0.36.0 From 7fd018d5eb3f44d28238c4c95dff4125c6b973a2 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 11 Dec 2023 13:18:58 +0100 Subject: [PATCH 8/8] Update Changelog --- CHANGELOG.md | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15d1a9c05..55f29a917 100644 --- a/CHANGELOG.md +++ b/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