Merge branch 'hotfix/0.35.1'

* hotfix/0.35.1: (24 commits) Update changelog Version 0.35.1 Add errors for healpix grid when ordering is not 'ring' Better error message when MatchingPartitioner is not initialised with Mesh/FunctionSpace to match atlas-grid-points Improve GridBuilder and atlas-grids output Fixup 03ef162 Use configurable KDTree geometry in PointCloud Configurable geometry in KDTree Fix warning about returning 0 instead of nullptr in ElementType.cc Fix warning about readability braces in ElementType.cc Fixes to MeshBuilder validate_mesh_vs_grid Fixup 0b1ed6b, fixing warnings of missing typename Add control to skip Gmsh-output of triangles with too large edge length ratio Add atlas::io support for std::vector<std::array> Better error message in atlas::io StdVectorAdaptor Allow constructor of atlas::io::ArrayShape with different integer types Support atlas_io with vector<std::int64_t> Make treat_failure_as_missing_value configurable Use search radius in FiniteElement when mesh defines metadata to do so ...
ecmwf · Oct 24, 2023 · 177f0ce · 177f0ce
2 parents 9b6fb72 + 7e6c7d2
commit 177f0ce
Show file tree

Hide file tree

Showing 36 changed files with 879 additions and 70 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,6 +7,23 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 
 ## [Unreleased]
 
+## [0.35.1] - 2023-24-10
+### 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
+- Configurable geometry in KDTree
+- Use configurable KDTree geometry in PointCloud
+- atlas-grid-points executable to list grid points
+- Allow constructor of atlas::io::ArrayShape with different integer types
+- Support atlas_io with vector<std::int64_t>
+
+### Fixed
+- Control FPE in StructuredColumns::checksum to avoid overflow and invalid in unimportant halo regions
+- Fixes to MeshBuilder validate_mesh_vs_grid
+
+### Changed
+- Use search radius in FiniteElement interpolation when mesh defines metadata to do so
+
 ## [0.35.0] - 2023-02-10
 ### Added
 - Add accessors for the GridBox class (MIR-632)
@@ -483,6 +500,7 @@ Fix StructuredInterpolation2D with retry for failed stencils
 ## 0.13.0 - 2018-02-16
 
 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
+[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
 [0.33.0]: https://github.com/ecmwf/atlas/compare/0.32.1...0.33.0

diff --git a/VERSION b/VERSION
@@ -1 +1 @@
-0.35.0
+0.35.1
diff --git a/atlas_io/src/atlas_io/CMakeLists.txt b/atlas_io/src/atlas_io/CMakeLists.txt
@@ -85,6 +85,7 @@ ecbuild_add_library( TARGET atlas_io
         types/array/ArrayReference.h
         types/array/adaptors/StdArrayAdaptor.h
         types/array/adaptors/StdVectorAdaptor.h
+        types/array/adaptors/StdVectorOfStdArrayAdaptor.h
         types/string.h
         types/scalar.h
         types/scalar.cc

diff --git a/atlas_io/src/atlas_io/detail/TypeTraits.h b/atlas_io/src/atlas_io/detail/TypeTraits.h
@@ -141,9 +141,11 @@ using enable_if_scalar_t = enable_if_t<std::is_scalar<T>::value && EnableBool>;
 
 template <typename T>
 constexpr bool is_array_datatype() {
-    return std::is_same<T, double>::value || std::is_same<T, float>::value ||
-                                std::is_same<T, int>::value || std::is_same<T, long>::value ||
-                                std::is_same<T, size_t>::value || std::is_same<T, std::byte>::value;
+    return std::is_same_v<T, double>        || std::is_same_v<T, float> ||
+           std::is_same_v<T, int>           || std::is_same_v<T, long>  ||
+           std::is_same_v<T, std::int32_t>  || std::is_same_v<T, std::int64_t> ||
+           std::is_same_v<T, std::uint64_t> || std::is_same_v<T, size_t> ||
+           std::is_same_v<T, std::byte>;
 }
 
 template <typename T>

diff --git a/atlas_io/src/atlas_io/types/array.h b/atlas_io/src/atlas_io/types/array.h
@@ -13,3 +13,4 @@
 #include "atlas_io/types/array/ArrayReference.h"
 #include "atlas_io/types/array/adaptors/StdArrayAdaptor.h"
 #include "atlas_io/types/array/adaptors/StdVectorAdaptor.h"
+#include "atlas_io/types/array/adaptors/StdVectorOfStdArrayAdaptor.h"
diff --git a/atlas_io/src/atlas_io/types/array/ArrayMetadata.h b/atlas_io/src/atlas_io/types/array/ArrayMetadata.h
@@ -12,6 +12,7 @@
 
 #include <cstdint>
 #include <string>
+#include <type_traits>
 
 #include "atlas_io/Metadata.h"
 #include "atlas_io/detail/DataType.h"
@@ -36,6 +37,32 @@ class ArrayShape : public std::vector<size_t> {
     ArrayShape(const std::array<idx_t, N>& list): Base(list.begin(), list.end()) {}
     template <typename idx_t>
     ArrayShape(const std::vector<idx_t>& list): Base(list.begin(), list.end()) {}
+    template <typename Int1, typename = std::enable_if_t<std::is_integral_v<Int1>>>
+    ArrayShape(Int1 i) {
+        resize(1);
+        operator[](0) = i;
+    }
+    template <typename Int1, typename Int2, typename = std::enable_if_t<std::is_integral_v<Int1> && std::is_integral_v<Int2>>>
+    ArrayShape(Int1 i, Int2 j) {
+        resize(2);
+        operator[](0) = i;
+        operator[](1) = j;
+    }
+    template <typename Int1, typename Int2, typename Int3, typename = std::enable_if_t<std::is_integral_v<Int1> && std::is_integral_v<Int2> && std::is_integral_v<Int3>>>
+    ArrayShape(Int1 i, Int2 j, Int3 k) {
+        resize(3);
+        operator[](0) = i;
+        operator[](1) = j;
+        operator[](2) = k;
+    }
+    template <typename Int1, typename Int2, typename Int3, typename Int4, typename = std::enable_if_t<std::is_integral_v<Int1> && std::is_integral_v<Int2> && std::is_integral_v<Int3> && std::is_integral_v<Int4>>>
+    ArrayShape(Int1 i, Int2 j, Int3 k, Int4 l) {
+        resize(4);
+        operator[](0) = i;
+        operator[](1) = j;
+        operator[](2) = k;
+        operator[](3) = l;
+    }
 };
 
 //---------------------------------------------------------------------------------------------------------------------

diff --git a/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h b/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h
@@ -37,7 +37,7 @@ void decode(const atlas::io::Metadata& m, const atlas::io::Data& encoded, std::v
     if (array.datatype().kind() != atlas::io::ArrayMetadata::DataType::kind<T>()) {
         std::stringstream err;
         err << "Could not decode " << m.json() << " into std::vector<" << atlas::io::demangle<T>() << ">. "
-            << "Incompatible datatype!";
+            << "Incompatible datatypes: " << array.datatype().str() << " and " << atlas::io::ArrayMetadata::DataType::str<T>() << ".";
         throw atlas::io::Exception(err.str(), Here());
     }
     const T* data = static_cast<const T*>(encoded.data());

diff --git a/atlas_io/src/atlas_io/types/array/adaptors/StdVectorOfStdArrayAdaptor.h b/atlas_io/src/atlas_io/types/array/adaptors/StdVectorOfStdArrayAdaptor.h
@@ -0,0 +1,70 @@
+/*
+ * (C) Copyright 2020 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 <array>
+
+#include "atlas_io/Data.h"
+#include "atlas_io/Exceptions.h"
+#include "atlas_io/Metadata.h"
+#include "atlas_io/types/array/ArrayMetadata.h"
+#include "atlas_io/types/array/ArrayReference.h"
+
+namespace std {
+
+//---------------------------------------------------------------------------------------------------------------------
+
+template <typename T, size_t N>
+void interprete(const std::vector<std::array<T, N>>& vector_of_array, atlas::io::ArrayReference& out) {
+    using atlas::io::ArrayReference;
+    out = ArrayReference{vector_of_array.front().data(), {vector_of_array.size(),N}};
+}
+
+//---------------------------------------------------------------------------------------------------------------------
+
+template <typename T, size_t N>
+void decode(const atlas::io::Metadata& m, const atlas::io::Data& encoded, std::vector<std::array<T, N>>& out) {
+    atlas::io::ArrayMetadata array(m);
+    if (array.datatype().kind() != atlas::io::ArrayMetadata::DataType::kind<T>()) {
+        std::stringstream err;
+        err << "Could not decode " << m.json() << " into std::vector<" << atlas::io::demangle<T>() << ">. "
+            << "Incompatible datatype!";
+        throw atlas::io::Exception(err.str(), Here());
+    }
+    if (array.rank() != 2) {
+        std::stringstream err;
+        err << "Could not decode " << m.json() << " into std::vector<std::array<" << atlas::io::demangle<T>() << "," << N << ">>. "
+            << "Incompatible rank!";
+        throw atlas::io::Exception(err.str(), Here());
+    }
+    if (array.shape(1) != N) {
+        std::stringstream err;
+        err << "Could not decode " << m.json() << " into std::vector<std::array<" << atlas::io::demangle<T>() << "," << N << ">>. "
+            << "Incompatible size!";
+        throw atlas::io::Exception(err.str(), Here());
+    }
+    const std::array<T,N>* data = static_cast<const std::array<T,N>*>(encoded.data());
+    // std::copy(data, data + array.shape(0), out.begin());
+    out.assign(data, data + array.shape(0));
+
+}
+
+//---------------------------------------------------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+}  // end namespace std
diff --git a/src/apps/CMakeLists.txt b/src/apps/CMakeLists.txt
@@ -24,6 +24,11 @@ ecbuild_add_executable(
   LIBS        atlas
   CONDITION   atlas_HAVE_ATLAS_GRID )
 
+ecbuild_add_executable(
+  TARGET      atlas-grid-points
+  SOURCES     atlas-grid-points.cc
+  LIBS        atlas )
+
 ecbuild_add_executable(
   TARGET      atlas-gaussian-latitudes
   SOURCES     atlas-gaussian-latitudes.cc

diff --git a/src/apps/atlas-grid-points.cc b/src/apps/atlas-grid-points.cc
@@ -0,0 +1,117 @@
+/*
+ * (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 <fstream>
+
+#include "atlas/grid.h"
+#include "atlas/library.h"
+#include "atlas/runtime/AtlasTool.h"
+#include "atlas/util/GridPointsJSONWriter.h"
+
+//------------------------------------------------------------------------------
+
+using namespace atlas;
+using eckit::PathName;
+
+//------------------------------------------------------------------------------
+
+class Program : public AtlasTool {
+    int execute(const Args& args) override;
+    std::string briefDescription() override { return "Write grid points to file"; }
+    std::string usage() override { return name() + " GRID [OPTION]... [--help]"; }
+    std::string longDescription() override {
+        return "    The 'GRID' argument can be either the name of a named grid, or the path to a"
+               " YAML configuration file that describes the grid.\n"
+               "Example values for grid names are: N80, F40, O24, L64x33, CS-ED-12. See the program "
+               "'atlas-grids' for a list of named grids.\n"
+               "\n";
+    }
+
+public:
+    Program(int argc, char** argv);
+
+private:
+};
+
+//-----------------------------------------------------------------------------
+
+Program::Program(int argc, char** argv): AtlasTool(argc, argv) {
+    add_option(new SimpleOption<std::string>("output.file", "Output file. If not specified, output is directed to stdout"));
+    add_option(new SimpleOption<std::string>("output.format", "Output format. If not specified: json"));
+    add_option(new SimpleOption<long>("field_base", "Base used for field output. Default=0"));
+    add_option(new SimpleOption<long>("index_base", "Base used for index input. Default=0"));
+    add_option(new SimpleOption<std::string>("index", 
+        "Select grid point indices (first_index=<index_base>, last_index = <index_base> + <size> - 1). "
+        "If not provided, all points are selected. "
+        "Format: comma separated list, where '-' can be used to represent a range. e.g. '[1-3,5,7-10]'."
+        "Square brackets are optional. white-spaces and newline characters are allowed as in a valid JSON array."));
+    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]"));
+    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"));
+    add_option(new SimpleOption<bool>("json.pretty", "Pretty printing of json output"));
+    add_option(new SimpleOption<bool>("verbose", "Output progress to stdout, default=false"));
+}
+
+//-----------------------------------------------------------------------------
+
+std::string get_arg(const AtlasTool::Args& args, const std::string& flag, const std::string& default_value = "") {
+    for (int i = 0; i < args.count() - 1; ++i) {
+        if (args(i) == flag) {
+            return args(i + 1);
+        }
+    }
+    if (not default_value.empty()) {
+        return default_value;
+    }
+    throw_Exception("Could not find argument for flag " + flag);
+}
+
+int Program::execute(const Args& args) {
+
+    Grid grid;
+    {
+        std::string key = args.count() ? args(0) : "";
+        if (!key.empty()) {
+            eckit::PathName path{key};
+            grid = path.exists() ? Grid(Grid::Spec{path}) : Grid(key);
+
+        }
+    }
+    if (!grid) {
+        Log::error() << "Grid not specified as positional argument" << std::endl;
+        return failed();
+    }
+
+    util::GridPointsJSONWriter writer{grid,args};
+    if (mpi::rank() == 0 ) {
+        std::string output_file;
+        if (args.get("output.file",output_file)) {
+            Log::info() << "Grid contains " << grid.size() << " points." << std::endl;
+            std::ofstream out(output_file);
+            writer.write(out, Log::info());
+            Log::info() << "File " << output_file << " written." << std::endl;
+        }
+        else {
+            writer.write(std::cout);
+        }
+    }
+    return success();
+}
+
+//------------------------------------------------------------------------------
+
+int main(int argc, char** argv) {
+    Program tool(argc, argv);
+    return tool.start();
+}
diff --git a/src/apps/atlas-grids.cc b/src/apps/atlas-grids.cc
@@ -146,8 +146,12 @@ int AtlasGrids::execute(const Args& args) {
         Log::info() << "usage: atlas-grids <grid> [OPTION]... [--help]\n" << std::endl;
         Log::info() << "\n";
         Log::info() << "Available grid types:" << std::endl;
-        for (auto b : grid::GridBuilder::typeRegistry()) {
-            Log::info() << "  -- " << b.second->type() << '\n';
+        std::set<std::string> grid_types;
+        for (const auto& [key, builder] : grid::GridBuilder::typeRegistry()) {
+            grid_types.insert(builder->type());
+        }
+        for (const auto& b : grid_types) {
+            Log::info() << "  -- " << b << '\n';
         }
         Log::info() << "\n";
         Log::info() << "Available named grids:" << std::endl;
@@ -159,13 +163,13 @@ int AtlasGrids::execute(const Args& args) {
             }
         }
 
-        for (auto b : grid::GridBuilder::nameRegistry()) {
+        for (const auto& [key, builder] : grid::GridBuilder::typeRegistry()) {
             int c = 0;
-            for (const auto& name : b.second->names()) {
+            for (const auto& name : builder->names()) {
                 if (c == 0) {
                     Log::info() << "  -- " << std::left << std::setw(maxlen + 8) << name;
-                    if (!b.second->type().empty()) {
-                        Log::info() << "type: " << b.second->type();
+                    if (!builder->type().empty()) {
+                        Log::info() << "type: " << builder->type();
                     }
                     Log::info() << std::endl;
                 }

diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt
@@ -799,6 +799,8 @@ util/DataType.h
 util/Earth.h
 util/Geometry.cc
 util/Geometry.h
+util/GridPointsJSONWriter.cc
+util/GridPointsJSONWriter.h
 util/KDTree.cc
 util/KDTree.h
 util/PolygonXY.cc

diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc
@@ -207,7 +207,7 @@ PointCloud::PointCloud(const Grid& grid, const grid::Partitioner& _partitioner,
             owned_lonlat.reserve(size_owned);
             owned_grid_idx.reserve(size_owned);
 
-            auto kdtree = util::IndexKDTree();
+            auto kdtree = util::IndexKDTree(config);
             {
                 ATLAS_TRACE("build kdtree");
                 kdtree.reserve(grid.size());

diff --git a/src/atlas/functionspace/detail/StructuredColumns.cc b/src/atlas/functionspace/detail/StructuredColumns.cc
@@ -27,6 +27,7 @@
 #include "atlas/grid/StructuredGrid.h"
 #include "atlas/grid/StructuredPartitionPolygon.h"
 #include "atlas/library/Library.h"
+#include "atlas/library/FloatingPointExceptions.h"
 #include "atlas/mesh/Mesh.h"
 #include "atlas/parallel/Checksum.h"
 #include "atlas/parallel/GatherScatter.h"
@@ -71,6 +72,8 @@ array::LocalView<T, 3> make_leveled_view(Field& field) {
 
 template <typename T>
 std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& field) {
+    bool disabled_fpe_overflow = library::disable_floating_point_exception(FE_OVERFLOW);
+    bool disabled_fpe_invalid  = library::disable_floating_point_exception(FE_INVALID);
     auto values = make_leveled_view<const T>(field);
     array::ArrayT<T> surface_field(values.shape(0), values.shape(2));
     auto surface     = array::make_view<T, 2>(surface_field);
@@ -83,6 +86,12 @@ std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& f
             }
         }
     }
+    if (disabled_fpe_overflow) {
+        library::enable_floating_point_exception(FE_OVERFLOW);
+    }
+    if (disabled_fpe_invalid) {
+        library::enable_floating_point_exception(FE_INVALID);
+    }
     return checksum.execute(surface.data(), surface_field.stride(0));
 }