From 6490e6a439b030d3b0e81d20a043c0be08020d82 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 5 Oct 2023 13:25:22 +0200 Subject: [PATCH 01/24] Update intel-oneapi install script --- tools/install-intel-oneapi.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/install-intel-oneapi.sh b/tools/install-intel-oneapi.sh index 78af1aec0..a954274d8 100755 --- a/tools/install-intel-oneapi.sh +++ b/tools/install-intel-oneapi.sh @@ -1,6 +1,6 @@ #!/bin/sh -KEY=GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB +KEY=GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB wget https://apt.repos.intel.com/intel-gpg-keys/$KEY sudo apt-key add $KEY rm $KEY From ad53b0f40c149fbeb1dcb15fabc9bf1596e260e3 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 6 Oct 2023 10:31:32 +0200 Subject: [PATCH 02/24] Control FPE in StructuredColumns::checksum to avoid overflow and invalid in unimportant halo regions --- src/atlas/functionspace/detail/StructuredColumns.cc | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/atlas/functionspace/detail/StructuredColumns.cc b/src/atlas/functionspace/detail/StructuredColumns.cc index 7ec33057b..653ea10ab 100644 --- 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 make_leveled_view(Field& field) { template 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(field); array::ArrayT surface_field(values.shape(0), values.shape(2)); auto surface = array::make_view(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)); } From eb696d217be3fe240ea342bf1f916985c25fa24d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Sat, 7 Oct 2023 00:06:27 +0200 Subject: [PATCH 03/24] Implement grid::Unstructured::footprint --- src/atlas/grid/detail/grid/Unstructured.cc | 5 +++++ src/atlas/grid/detail/grid/Unstructured.h | 2 ++ 2 files changed, 7 insertions(+) diff --git a/src/atlas/grid/detail/grid/Unstructured.cc b/src/atlas/grid/detail/grid/Unstructured.cc index 258436b6e..b18e8de09 100644 --- a/src/atlas/grid/detail/grid/Unstructured.cc +++ b/src/atlas/grid/detail/grid/Unstructured.cc @@ -182,6 +182,11 @@ void Unstructured::hash(eckit::Hash& h) const { projection().hash(h); } +size_t Unstructured::footprint() const { + return sizeof(PointXY) * points_->size(); +} + + RectangularLonLatDomain Unstructured::lonlatBoundingBox() const { return projection_ ? projection_.lonlatBoundingBox(domain_) : domain_; } diff --git a/src/atlas/grid/detail/grid/Unstructured.h b/src/atlas/grid/detail/grid/Unstructured.h index 9c830d22e..e99942f69 100644 --- a/src/atlas/grid/detail/grid/Unstructured.h +++ b/src/atlas/grid/detail/grid/Unstructured.h @@ -196,6 +196,8 @@ class Unstructured : public Grid { Config meshgenerator() const override; Config partitioner() const override; + size_t footprint() const override; + private: // methods virtual void print(std::ostream&) const override; From 0ad41fe786c739229a16297bc6601df0600f5859 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 9 Oct 2023 11:37:07 +0200 Subject: [PATCH 04/24] Don't output with output::Gmsh the triangle elements with wrong orientation when coordinates are "lonlat" --- src/atlas/output/detail/GmshIO.cc | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/atlas/output/detail/GmshIO.cc b/src/atlas/output/detail/GmshIO.cc index 8a3979b16..9a5dd3736 100644 --- a/src/atlas/output/detail/GmshIO.cc +++ b/src/atlas/output/detail/GmshIO.cc @@ -913,6 +913,7 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { bool coords_is_idx = coords_field.datatype().kind() == array::make_datatype().kind(); auto coords = array::make_view(coords_is_idx ? dummy_double : coords_field.array()); auto coords_idx = array::make_view(coords_is_idx ? coords_field.array() : dummy_idx); + bool coords_is_lonlat = coords_field.name() == "lonlat"; auto glb_idx = array::make_view(nodes.global_index()); @@ -989,6 +990,7 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { nb_elements += hybrid->size(); const auto hybrid_halo = array::make_view(hybrid->halo()); const auto hybrid_flags = array::make_view(hybrid->flags()); + const auto& node_connectivity = hybrid->node_connectivity(); auto include = [&](idx_t e) { auto topology = Topology::view(hybrid_flags(e)); @@ -1013,6 +1015,21 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { if (topology.check(Topology::INVALID)) { return false; } + + if (coords_is_lonlat) { + if (node_connectivity.cols(e) == 3) { + auto x0 = [&]() { return coords(node_connectivity(e,0),LON); }; + auto x1 = [&]() { return coords(node_connectivity(e,1),LON); }; + auto x2 = [&]() { return coords(node_connectivity(e,2),LON); }; + auto y0 = [&]() { return coords(node_connectivity(e,0),LAT); }; + auto y1 = [&]() { return coords(node_connectivity(e,1),LAT); }; + auto y2 = [&]() { return coords(node_connectivity(e,2),LAT); }; + auto triangle_area = (x0() * (y1() - y2()) + x1() * (y2() - y0()) + x2() * (y0() - y1())) * 0.5; + if (triangle_area <= 0 ) { + return false; + } + } + } return true; }; @@ -1081,6 +1098,20 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { if (topology.check(Topology::INVALID)) { return false; } + if (coords_is_lonlat) { + if (nb_nodes == 3) { + auto x0 = [&]() { return coords(node_connectivity(e,0),LON); }; + auto x1 = [&]() { return coords(node_connectivity(e,1),LON); }; + auto x2 = [&]() { return coords(node_connectivity(e,2),LON); }; + auto y0 = [&]() { return coords(node_connectivity(e,0),LAT); }; + auto y1 = [&]() { return coords(node_connectivity(e,1),LAT); }; + auto y2 = [&]() { return coords(node_connectivity(e,2),LAT); }; + auto triangle_area = (x0() * (y1() - y2()) + x1() * (y2() - y0()) + x2() * (y0() - y1())) * 0.5; + if (triangle_area <= 0 ) { + return false; + } + } + } return true; }; if (binary) { From c728ae5923ba873c18cba2592770b0137fa310e7 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Sat, 7 Oct 2023 00:08:34 +0200 Subject: [PATCH 05/24] Use search radius in FiniteElement when mesh defines metadata to do so --- .../method/unstructured/FiniteElement.cc | 37 ++++++++++++++----- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/src/atlas/interpolation/method/unstructured/FiniteElement.cc b/src/atlas/interpolation/method/unstructured/FiniteElement.cc index 74c3e63c4..3e41003da 100644 --- a/src/atlas/interpolation/method/unstructured/FiniteElement.cc +++ b/src/atlas/interpolation/method/unstructured/FiniteElement.cc @@ -250,6 +250,13 @@ void FiniteElement::setup(const FunctionSpace& source) { const idx_t maxNbElemsToTry = std::max(8, idx_t(Nelements * max_fraction_elems_to_try_)); idx_t max_neighbours = 0; + double search_radius = 0.; + if (meshSource.metadata().has("cell_maximum_diagonal_on_unit_sphere")) { + search_radius = Geometry("Earth").radius() * meshSource.metadata().getDouble("cell_maximum_diagonal_on_unit_sphere"); + ASSERT(search_radius > 0.); + Log::debug() << "k-d tree: search radius = " << search_radius/1000. << " km" << std::endl; + } + std::vector failures; ATLAS_TRACE_SCOPE("Computing interpolation matrix") { @@ -265,17 +272,29 @@ void FiniteElement::setup(const FunctionSpace& source) { bool success = false; std::ostringstream failures_log; - while (!success && kpts <= maxNbElemsToTry) { - max_neighbours = std::max(kpts, max_neighbours); + if (search_radius != 0.) { + ElemIndex3::NodeList cs = eTree->findInSphere(p,search_radius); + if (cs.size()) { + Triplets triplets = projectPointToElements(ip, cs, failures_log); - ElemIndex3::NodeList cs = eTree->kNearestNeighbours(p, kpts); - Triplets triplets = projectPointToElements(ip, cs, failures_log); - - if (triplets.size()) { - std::copy(triplets.begin(), triplets.end(), std::back_inserter(weights_triplets)); - success = true; + if (triplets.size()) { + std::copy(triplets.begin(), triplets.end(), std::back_inserter(weights_triplets)); + success = true; + } + } + } + else { + while (!success && kpts <= maxNbElemsToTry) { + max_neighbours = std::max(kpts, max_neighbours); + ElemIndex3::NodeList cs = eTree->kNearestNeighbours(p, kpts); + Triplets triplets = projectPointToElements(ip, cs, failures_log); + + if (triplets.size()) { + std::copy(triplets.begin(), triplets.end(), std::back_inserter(weights_triplets)); + success = true; + } + kpts *= 2; } - kpts *= 2; } if (!success) { From dd3987547a6aefb7f3c29d1ff937d64e56331d8f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Sat, 7 Oct 2023 00:09:46 +0200 Subject: [PATCH 06/24] Make treat_failure_as_missing_value configurable --- src/atlas/interpolation/method/unstructured/FiniteElement.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/atlas/interpolation/method/unstructured/FiniteElement.h b/src/atlas/interpolation/method/unstructured/FiniteElement.h index 0b7d99458..8fe5bc08e 100644 --- a/src/atlas/interpolation/method/unstructured/FiniteElement.h +++ b/src/atlas/interpolation/method/unstructured/FiniteElement.h @@ -30,6 +30,7 @@ class FiniteElement : public Method { public: FiniteElement(const Config& config): Method(config) { config.get("max_fraction_elems_to_try", max_fraction_elems_to_try_); + config.get("treat_failure_as_missing_value", treat_failure_as_missing_value_); } virtual ~FiniteElement() override {} From 24b8386319677693410eab9aed82c0eb65307e13 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Sat, 7 Oct 2023 00:25:56 +0200 Subject: [PATCH 07/24] Support atlas_io with vector --- atlas_io/src/atlas_io/detail/TypeTraits.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/atlas_io/src/atlas_io/detail/TypeTraits.h b/atlas_io/src/atlas_io/detail/TypeTraits.h index 2ecf759d3..18c1e1c3e 100644 --- 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::value && EnableBool>; template constexpr bool is_array_datatype() { - return std::is_same::value || std::is_same::value || - std::is_same::value || std::is_same::value || - std::is_same::value || std::is_same::value; + return std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v || + std::is_same_v; } template From 0b1ed6b060947da8e6f7df6b104a1d30c85bf0fc Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 9 Oct 2023 11:41:54 +0200 Subject: [PATCH 08/24] Allow constructor of atlas::io::ArrayShape with different integer types --- .../src/atlas_io/types/array/ArrayMetadata.h | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/atlas_io/src/atlas_io/types/array/ArrayMetadata.h b/atlas_io/src/atlas_io/types/array/ArrayMetadata.h index 4711afb44..19aacdcd3 100644 --- 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 #include +#include #include "atlas_io/Metadata.h" #include "atlas_io/detail/DataType.h" @@ -36,6 +37,32 @@ class ArrayShape : public std::vector { ArrayShape(const std::array& list): Base(list.begin(), list.end()) {} template ArrayShape(const std::vector& list): Base(list.begin(), list.end()) {} + template >> + ArrayShape(Int1 i) { + resize(1); + operator[](0) = i; + } + template && std::is_integral_v>> + ArrayShape(Int1 i, Int2 j) { + resize(2); + operator[](0) = i; + operator[](1) = j; + } + template && std::is_integral_v && std::is_integral_v>> + ArrayShape(Int1 i, Int2 j, Int3 k) { + resize(3); + operator[](0) = i; + operator[](1) = j; + operator[](2) = k; + } + template && std::is_integral_v && std::is_integral_v && std::is_integral_v>> + ArrayShape(Int1 i, Int2 j, Int3 k, Int4 l) { + resize(4); + operator[](0) = i; + operator[](1) = j; + operator[](2) = k; + operator[](3) = l; + } }; //--------------------------------------------------------------------------------------------------------------------- From 128ab3e5711046a37848fab76aa56e59541050b2 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 9 Oct 2023 11:43:04 +0200 Subject: [PATCH 09/24] Better error message in atlas::io StdVectorAdaptor --- atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h b/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h index b2be76a32..1b2859e26 100644 --- 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()) { std::stringstream err; err << "Could not decode " << m.json() << " into std::vector<" << atlas::io::demangle() << ">. " - << "Incompatible datatype!"; + << "Incompatible datatypes: " << array.datatype().str() << " and " << atlas::io::ArrayMetadata::DataType::str() << "."; throw atlas::io::Exception(err.str(), Here()); } const T* data = static_cast(encoded.data()); From a14961fc3d81eed8f91546706bb0c59812a26cfc Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 9 Oct 2023 11:47:47 +0200 Subject: [PATCH 10/24] Add atlas::io support for std::vector --- atlas_io/src/atlas_io/CMakeLists.txt | 1 + atlas_io/src/atlas_io/types/array.h | 1 + .../adaptors/StdVectorOfStdArrayAdaptor.h | 70 +++++++++++++++++++ 3 files changed, 72 insertions(+) create mode 100644 atlas_io/src/atlas_io/types/array/adaptors/StdVectorOfStdArrayAdaptor.h diff --git a/atlas_io/src/atlas_io/CMakeLists.txt b/atlas_io/src/atlas_io/CMakeLists.txt index 395d827d2..877d61ead 100644 --- 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/types/array.h b/atlas_io/src/atlas_io/types/array.h index d2129e674..f205d91bb 100644 --- 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/adaptors/StdVectorOfStdArrayAdaptor.h b/atlas_io/src/atlas_io/types/array/adaptors/StdVectorOfStdArrayAdaptor.h new file mode 100644 index 000000000..1949b6c2b --- /dev/null +++ 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 +#include + +#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 +void interprete(const std::vector>& vector_of_array, atlas::io::ArrayReference& out) { + using atlas::io::ArrayReference; + out = ArrayReference{vector_of_array.front().data(), {vector_of_array.size(),N}}; +} + +//--------------------------------------------------------------------------------------------------------------------- + +template +void decode(const atlas::io::Metadata& m, const atlas::io::Data& encoded, std::vector>& out) { + atlas::io::ArrayMetadata array(m); + if (array.datatype().kind() != atlas::io::ArrayMetadata::DataType::kind()) { + std::stringstream err; + err << "Could not decode " << m.json() << " into std::vector<" << atlas::io::demangle() << ">. " + << "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() << "," << 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() << "," << N << ">>. " + << "Incompatible size!"; + throw atlas::io::Exception(err.str(), Here()); + } + const std::array* data = static_cast*>(encoded.data()); + // std::copy(data, data + array.shape(0), out.begin()); + out.assign(data, data + array.shape(0)); + +} + +//--------------------------------------------------------------------------------------------------------------------- + + + + + + + + +} // end namespace std From d25be1b71ac19bd26d3f44f44a8f4d147c8f5f5d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 12:44:33 +0200 Subject: [PATCH 11/24] Add control to skip Gmsh-output of triangles with too large edge length ratio --- src/atlas/output/detail/GmshIO.cc | 37 +++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/atlas/output/detail/GmshIO.cc b/src/atlas/output/detail/GmshIO.cc index 9a5dd3736..305051de4 100644 --- a/src/atlas/output/detail/GmshIO.cc +++ b/src/atlas/output/detail/GmshIO.cc @@ -17,6 +17,7 @@ #include #include "eckit/filesystem/PathName.h" +#include "eckit/config/Resource.h" #include "atlas/array.h" #include "atlas/array/ArrayView.h" @@ -915,6 +916,8 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { auto coords_idx = array::make_view(coords_is_idx ? coords_field.array() : dummy_idx); bool coords_is_lonlat = coords_field.name() == "lonlat"; + double filter_edge_ratio = eckit::Resource("$ATLAS_GMSH_FILTER_EDGE_RATIO",0.); + auto glb_idx = array::make_view(nodes.global_index()); const idx_t surfdim = nodes.field(nodes_field).shape(1); // nb of variables in coords @@ -1028,6 +1031,23 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { if (triangle_area <= 0 ) { return false; } + // edge length comparison + if (filter_edge_ratio > 0.) { + auto dx10 = (x1()-x0()); + auto dx21 = (x2()-x1()); + auto dx02 = (x0()-x2()); + auto dy10 = (y1()-y0()); + auto dy21 = (y2()-y1()); + auto dy02 = (y0()-y2()); + auto d10 = dx10*dx10 + dy10*dy10; + auto d21 = dx21*dx21 + dy21*dy21; + auto d02 = dx02*dx02 + dy02*dy02; + auto dmin = std::min(d10, std::min(d21, d02)); + auto dmax = std::max(d10, std::max(d21, d02)); + if (dmax > filter_edge_ratio * dmin) { + return false; + } + } } } return true; @@ -1110,6 +1130,23 @@ void GmshIO::write(const Mesh& mesh, const PathName& file_path) const { if (triangle_area <= 0 ) { return false; } + // edge length comparison + if (filter_edge_ratio > 0.) { + auto dx10 = (x1()-x0()); + auto dx21 = (x2()-x1()); + auto dx02 = (x0()-x2()); + auto dy10 = (y1()-y0()); + auto dy21 = (y2()-y1()); + auto dy02 = (y0()-y2()); + auto d10 = dx10*dx10 + dy10*dy10; + auto d21 = dx21*dx21 + dy21*dy21; + auto d02 = dx02*dx02 + dy02*dy02; + auto dmin = std::min(d10, std::min(d21, d02)); + auto dmax = std::max(d10, std::max(d21, d02)); + if (dmax > filter_edge_ratio * dmin) { + return false; + } + } } } return true; From b4108aafa93fd0b330052e9e45b759bfbd95215f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 11:01:51 +0000 Subject: [PATCH 12/24] Fixup 0b1ed6b06, fixing warnings of missing typename --- atlas_io/src/atlas_io/types/array/ArrayMetadata.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/atlas_io/src/atlas_io/types/array/ArrayMetadata.h b/atlas_io/src/atlas_io/types/array/ArrayMetadata.h index 19aacdcd3..bd0015285 100644 --- a/atlas_io/src/atlas_io/types/array/ArrayMetadata.h +++ b/atlas_io/src/atlas_io/types/array/ArrayMetadata.h @@ -37,25 +37,25 @@ class ArrayShape : public std::vector { ArrayShape(const std::array& list): Base(list.begin(), list.end()) {} template ArrayShape(const std::vector& list): Base(list.begin(), list.end()) {} - template >> + template >> ArrayShape(Int1 i) { resize(1); operator[](0) = i; } - template && std::is_integral_v>> + template && std::is_integral_v>> ArrayShape(Int1 i, Int2 j) { resize(2); operator[](0) = i; operator[](1) = j; } - template && std::is_integral_v && std::is_integral_v>> + template && std::is_integral_v && std::is_integral_v>> ArrayShape(Int1 i, Int2 j, Int3 k) { resize(3); operator[](0) = i; operator[](1) = j; operator[](2) = k; } - template && std::is_integral_v && std::is_integral_v && std::is_integral_v>> + template && std::is_integral_v && std::is_integral_v && std::is_integral_v>> ArrayShape(Int1 i, Int2 j, Int3 k, Int4 l) { resize(4); operator[](0) = i; From 03ef162cba532522e546db166f2eeb426c97805b Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 11:35:16 +0000 Subject: [PATCH 13/24] Fixes to MeshBuilder validate_mesh_vs_grid - comparison function was multiplying instead of adding 1.e-16 - scoped index i was shadowing loop index i --- src/atlas/mesh/MeshBuilder.cc | 46 +++++++++++++++++------------------ 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/atlas/mesh/MeshBuilder.cc b/src/atlas/mesh/MeshBuilder.cc index 72a6c6f41..735eb2e2d 100644 --- a/src/atlas/mesh/MeshBuilder.cc +++ b/src/atlas/mesh/MeshBuilder.cc @@ -31,11 +31,11 @@ atlas::UnstructuredGrid assemble_unstructured_grid(size_t nb_nodes, const double const size_t nb_owned_nodes = std::count(ghosts, ghosts + nb_nodes, 0); std::vector owned_lonlats(2 * nb_owned_nodes); int counter = 0; - for (size_t i = 0; i < nb_nodes; ++i) { - if (ghosts[i] == 0) { - owned_lonlats[counter] = lons[i]; + for (size_t n = 0; n < nb_nodes; ++n) { + if (ghosts[n] == 0) { + owned_lonlats[counter] = lons[n]; counter++; - owned_lonlats[counter] = lats[i]; + owned_lonlats[counter] = lats[n]; counter++; } } @@ -51,8 +51,8 @@ atlas::UnstructuredGrid assemble_unstructured_grid(size_t nb_nodes, const double global_lonlats = std::move(buffer.buffer); std::vector points(nb_nodes_global); - for (size_t i = 0; i < nb_nodes_global; ++i) { - points[i] = atlas::PointXY({global_lonlats[2*i], global_lonlats[2*i + 1]}); + for (size_t n = 0; n < nb_nodes_global; ++n) { + points[n] = atlas::PointXY({global_lonlats[2*n], global_lonlats[2*n + 1]}); } return atlas::UnstructuredGrid(new std::vector(points.begin(), points.end())); @@ -73,47 +73,47 @@ void validate_grid_vs_mesh(const atlas::Grid& grid, size_t nb_nodes, const doubl const size_t nb_owned_nodes = std::count(ghosts, ghosts + nb_nodes, 0); size_t nb_nodes_global = 0; comm.allReduce(nb_owned_nodes, nb_nodes_global, eckit::mpi::sum()); - for (size_t i = 0; i < nb_nodes; ++i) { - if (ghosts[i] == 0) { + for (size_t n = 0; n < nb_nodes; ++n) { + if (ghosts[n] == 0) { // Check global_indices is consistent with a 1-based contiguous index over nodes - ATLAS_ASSERT(global_indices[i] >= 1); - ATLAS_ASSERT(global_indices[i] <= nb_nodes_global); + ATLAS_ASSERT(global_indices[n] >= 1); + ATLAS_ASSERT(global_indices[n] <= nb_nodes_global); } } double lonlat[2]; const auto equal_within_roundoff = [](const double a, const double b) -> bool { - return fabs(a - b) <= 360.0 * 1.0e-16; + return std::abs(a - b) <= 360.0 + 1.0e-16; }; // Check lonlats for each supported grid type if (grid.type() == "unstructured") { const atlas::UnstructuredGrid ugrid(grid); - for (size_t i = 0; i < nb_nodes; ++i) { - if (ghosts[i] == 0) { - ugrid.lonlat(global_indices[i] - 1, lonlat); - if (!equal_within_roundoff(lonlat[0], lons[i]) || !equal_within_roundoff(lonlat[1], lats[i])) { + for (size_t n = 0; n < nb_nodes; ++n) { + if (ghosts[n] == 0) { + ugrid.lonlat(global_indices[n] - 1, lonlat); + if (!equal_within_roundoff(lonlat[0], lons[n]) || !equal_within_roundoff(lonlat[1], lats[n])) { throw_Exception("In MeshBuilder: UnstructuredGrid from config does not match mesh coordinates", Here()); } } } } else if (grid.type() == "structured") { const atlas::StructuredGrid sgrid(grid); - for (size_t i = 0; i < nb_nodes; ++i) { - if (ghosts[i] == 0) { + for (size_t n = 0; n < nb_nodes; ++n) { + if (ghosts[n] == 0) { idx_t i, j; - sgrid.index2ij(global_indices[i] - 1, i, j); + sgrid.index2ij(global_indices[n] - 1, i, j); sgrid.lonlat(i, j, lonlat); - if (!equal_within_roundoff(lonlat[0], lons[i]) || !equal_within_roundoff(lonlat[1], lats[i])) { + if (!equal_within_roundoff(lonlat[0], lons[n]) || !equal_within_roundoff(lonlat[1], lats[n])) { throw_Exception("In MeshBuilder: StructuredGrid from config does not match mesh coordinates", Here()); } } } } else { - for (size_t i = 0; i < nb_nodes; ++i) { - if (ghosts[i] == 0) { - auto point = *(grid.lonlat().begin() + static_cast(global_indices[i] - 1)); - if (!equal_within_roundoff(point.lon(), lons[i]) || !equal_within_roundoff(point.lat(), lats[i])) { + for (size_t n = 0; n < nb_nodes; ++n) { + if (ghosts[n] == 0) { + auto point = *(grid.lonlat().begin() + static_cast(global_indices[n] - 1)); + if (!equal_within_roundoff(point.lon(), lons[n]) || !equal_within_roundoff(point.lat(), lats[n])) { throw_Exception("In MeshBuilder: Grid from config does not match mesh coordinates", Here()); } } From 4840b9cc6c91068be6f6697078a822aa80be3f6f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 14:21:53 +0200 Subject: [PATCH 14/24] Fix warning about readability braces in ElementType.cc --- src/atlas/mesh/ElementType.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/atlas/mesh/ElementType.cc b/src/atlas/mesh/ElementType.cc index 60c033d3d..b027f80fa 100644 --- a/src/atlas/mesh/ElementType.cc +++ b/src/atlas/mesh/ElementType.cc @@ -25,15 +25,15 @@ ElementType* ElementType::create( const std::string& name ) { // currently naive implementation. This has to be replaced // with self-registration and factory mechanism - if( name == "Triangle" ) + if( name == "Triangle" ) { return new elementtypes::Triangle(); - else if( name == "Quadrilateral") + } else if( name == "Quadrilateral") { return new elementtypes::Quadrilateral(); - else if( name == "Line" ) + } else if( name == "Line" ) { return new elementtypes::Line(); - else if( name == "Pentagon") + } else if( name == "Pentagon") { return new elementtypes::Pentagon(); - else { + } else { throw_NotImplemented("Element type ["+name+"] does not exist", Here()); } return 0; From 0738038af9d58c8affcaffd82940a6fa09f06053 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 15:59:20 +0200 Subject: [PATCH 15/24] Fix warning about returning 0 instead of nullptr in ElementType.cc --- src/atlas/mesh/ElementType.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/mesh/ElementType.cc b/src/atlas/mesh/ElementType.cc index b027f80fa..4865b52d9 100644 --- a/src/atlas/mesh/ElementType.cc +++ b/src/atlas/mesh/ElementType.cc @@ -36,7 +36,7 @@ ElementType* ElementType::create( const std::string& name ) } else { throw_NotImplemented("Element type ["+name+"] does not exist", Here()); } - return 0; + return nullptr; } ElementType::ElementType() = default; From a89b7871975d8ba4432a4557e5009830289ef516 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 16:52:43 +0200 Subject: [PATCH 16/24] Configurable geometry in KDTree --- src/atlas/util/KDTree.h | 6 ++++++ src/tests/util/test_kdtree.cc | 22 ++++++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/src/atlas/util/KDTree.h b/src/atlas/util/KDTree.h index 898f0d813..b31936304 100644 --- a/src/atlas/util/KDTree.h +++ b/src/atlas/util/KDTree.h @@ -10,6 +10,8 @@ #pragma once +#include "eckit/config/Configuration.h" + #include "atlas/util/Geometry.h" #include "atlas/util/ObjectHandle.h" #include "atlas/util/detail/KDTree.h" @@ -71,6 +73,10 @@ class KDTree : public ObjectHandle> { /// @brief Construct an empty kd-tree with custom geometry KDTree(const Geometry& geometry): Handle(new detail::KDTreeMemory(geometry)) {} + /// @brief Construct an empty kd-tree with custom geometry + KDTree(const eckit::Configuration& config): + KDTree(Geometry(config.getString("geometry","Earth"))) {} + /// @brief Construct a shared kd-tree with default geometry (Earth) template KDTree(const std::shared_ptr& kdtree): Handle(new detail::KDTree_eckit(kdtree)) {} diff --git a/src/tests/util/test_kdtree.cc b/src/tests/util/test_kdtree.cc index 782a0f7b6..073324a84 100644 --- a/src/tests/util/test_kdtree.cc +++ b/src/tests/util/test_kdtree.cc @@ -180,6 +180,7 @@ CASE("test kdtree") { EXPECT_EQ(neighbours, expected_neighbours); } + CASE("test assertion") { auto grid = Grid{"O32"}; @@ -312,6 +313,27 @@ CASE("test IndexKDTree 2D vs 3D") { // Note that the expected values are different whether 2D search or 3D search is used } + +CASE("test kdtree with configured geometry") { + auto grid = Grid{"O32"}; + + IndexKDTree search_unit(util::Config("geometry","UnitSphere")); + IndexKDTree search_earth(util::Config("geometry","Earth")); + + search_unit.build(grid.lonlat(),PayloadGenerator(grid.size())); + search_earth.build(grid.lonlat(),PayloadGenerator(grid.size())); + + double km_unit = 1000. / util::Earth::radius(); + double km_earth = 1000.; + auto neighbours_unit = search_unit. closestPointsWithinRadius(PointLonLat{180., 45.}, 500 * km_unit ).payloads(); + auto neighbours_earth = search_earth.closestPointsWithinRadius(PointLonLat{180., 45.}, 500 * km_earth).payloads(); + + auto expected_neighbours = std::vector{760, 842, 759, 761, 841, 843, 682}; + + EXPECT_EQ(neighbours_unit, expected_neighbours); + EXPECT_EQ(neighbours_earth, expected_neighbours); +} + //------------------------------------------------------------------------------------------------ } // namespace test From 1fa9822859bc7d355e20180f15609a52785e6176 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 13 Oct 2023 17:25:22 +0200 Subject: [PATCH 17/24] Use configurable KDTree geometry in PointCloud --- src/atlas/functionspace/PointCloud.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc index 4083d60d1..61388a2f5 100644 --- 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()); From b8a595ad9ead53e75a0feee6d1d0175c67074ee8 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 16 Oct 2023 11:17:20 +0200 Subject: [PATCH 18/24] Fixup 03ef162 --- src/atlas/mesh/MeshBuilder.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/mesh/MeshBuilder.cc b/src/atlas/mesh/MeshBuilder.cc index 735eb2e2d..32020d997 100644 --- a/src/atlas/mesh/MeshBuilder.cc +++ b/src/atlas/mesh/MeshBuilder.cc @@ -83,7 +83,7 @@ void validate_grid_vs_mesh(const atlas::Grid& grid, size_t nb_nodes, const doubl double lonlat[2]; const auto equal_within_roundoff = [](const double a, const double b) -> bool { - return std::abs(a - b) <= 360.0 + 1.0e-16; + return std::abs(a - b) <= 360.0 * 1.0e-16; }; // Check lonlats for each supported grid type From e38becb86cffb9bd86a842ab14eac6ee7814973f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Oct 2023 09:22:04 +0200 Subject: [PATCH 19/24] Improve GridBuilder and atlas-grids output Adding GridBuilder::registerNamedGrid, which allows to add grids in a plugin programatically --- src/apps/atlas-grids.cc | 16 ++++++++++------ src/atlas/grid/detail/grid/Grid.cc | 12 ++++++------ src/atlas/grid/detail/grid/GridBuilder.cc | 10 ++++++++++ src/atlas/grid/detail/grid/GridBuilder.h | 4 +++- 4 files changed, 29 insertions(+), 13 deletions(-) diff --git a/src/apps/atlas-grids.cc b/src/apps/atlas-grids.cc index 953795985..bca70e95d 100644 --- 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 [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 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/grid/detail/grid/Grid.cc b/src/atlas/grid/detail/grid/Grid.cc index 467076d39..e9750b3d7 100644 --- a/src/atlas/grid/detail/grid/Grid.cc +++ b/src/atlas/grid/detail/grid/Grid.cc @@ -71,9 +71,8 @@ const Grid* Grid::create(const std::string& name) { } const Grid* Grid::create(const std::string& name, const Grid::Config& config) { - const GridBuilder::Registry& registry = GridBuilder::nameRegistry(); - for (GridBuilder::Registry::const_iterator it = registry.begin(); it != registry.end(); ++it) { - const Grid* grid = it->second->create(name, config); + for (const auto& [key, builder]: GridBuilder::nameRegistry()) { + const Grid* grid = builder->create(name, config); if (grid) { return grid; } @@ -83,11 +82,12 @@ const Grid* Grid::create(const std::string& name, const Grid::Config& config) { std::ostringstream log; log << "Could not construct Grid from the name \"" << name << "\"\n"; log << "Accepted names are: \n"; - for (GridBuilder::Registry::const_iterator it = registry.begin(); it != registry.end(); ++it) { - log << " - " << *it->second << "\n"; + for (const auto& [key, grid_builder]: GridBuilder::typeRegistry()) { + for( auto& grid_name: grid_builder->names()) { + log << " - " << grid_name << "\n"; + } } throw_Exception(log.str()); - // return GridBuilder::createNamed(name); } const Grid* Grid::create(const Grid& grid, const Domain& domain) { diff --git a/src/atlas/grid/detail/grid/GridBuilder.cc b/src/atlas/grid/detail/grid/GridBuilder.cc index 06166506d..acb8276ed 100644 --- a/src/atlas/grid/detail/grid/GridBuilder.cc +++ b/src/atlas/grid/detail/grid/GridBuilder.cc @@ -179,6 +179,16 @@ GridBuilder::~GridBuilder() { } } +void GridBuilder::registerNamedGrid(const std::string& name) { + pthread_once(&once, init); + eckit::AutoLock lock(local_mutex); + std::string regex = "^"+name+"$"; + ATLAS_ASSERT(named_grids->find(regex) == named_grids->end()); + (*named_grids)[regex] = this; + names_.emplace_back(regex); + pretty_names_.emplace_back(name); +} + const Grid::Implementation* GridBuilder::create(const Grid::Config& config) const { //eckit::Factory& fact = eckit::Factory::instance(); diff --git a/src/atlas/grid/detail/grid/GridBuilder.h b/src/atlas/grid/detail/grid/GridBuilder.h index 5e8454ce8..8b7d754aa 100644 --- a/src/atlas/grid/detail/grid/GridBuilder.h +++ b/src/atlas/grid/detail/grid/GridBuilder.h @@ -39,10 +39,12 @@ class GridBuilder { virtual const Grid::Implementation* create(const std::string&, const Grid::Config& = Grid::Config()) const = 0; - const std::string& type() const; + virtual const std::string& type() const; const std::vector& regexes() const; const std::vector& names() const; + void registerNamedGrid(const std::string& name); + protected: bool match(const std::string& string, std::vector& matches, int& id) const; From 22f02c0284b905d28c967c3a59ca06f3fbcd4782 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 17 Apr 2023 13:38:37 -0600 Subject: [PATCH 20/24] atlas-grid-points --- src/apps/CMakeLists.txt | 5 + src/apps/atlas-grid-points.cc | 117 +++++++++ src/atlas/CMakeLists.txt | 2 + src/atlas/util/GridPointsJSONWriter.cc | 339 +++++++++++++++++++++++++ src/atlas/util/GridPointsJSONWriter.h | 48 ++++ 5 files changed, 511 insertions(+) create mode 100644 src/apps/atlas-grid-points.cc create mode 100644 src/atlas/util/GridPointsJSONWriter.cc create mode 100644 src/atlas/util/GridPointsJSONWriter.h diff --git a/src/apps/CMakeLists.txt b/src/apps/CMakeLists.txt index 2e8eb08ac..163be727f 100644 --- 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 new file mode 100644 index 000000000..86cf4c674 --- /dev/null +++ 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 + +#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("output.file", "Output file. If not specified, output is directed to stdout")); + add_option(new SimpleOption("output.format", "Output format. If not specified: json")); + add_option(new SimpleOption("field_base", "Base used for field output. Default=0")); + add_option(new SimpleOption("index_base", "Base used for index input. Default=0")); + add_option(new SimpleOption("index", + "Select grid point indices (first_index=, last_index = + - 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("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]")); + 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")); + add_option(new SimpleOption("json.pretty", "Pretty printing of json output")); + add_option(new SimpleOption("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/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index dad416859..2b92a933b 100644 --- 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/util/GridPointsJSONWriter.cc b/src/atlas/util/GridPointsJSONWriter.cc new file mode 100644 index 000000000..e7066ad77 --- /dev/null +++ b/src/atlas/util/GridPointsJSONWriter.cc @@ -0,0 +1,339 @@ +/* + * (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. + */ + +#include "GridPointsJSONWriter.h" + +#include +#include +#include +#include + +#include "eckit/utils/Tokenizer.h" + +namespace atlas { +namespace util { + +//------------------------------------------------------------------------------ + +namespace { +std::vector points_from_list(const std::string& list, long base) { + std::vector points; + if (list.empty()) { + return points; + } + if (list[0] == '[' && list[list.size()-1] == ']') { + return points_from_list(std::string(list.begin()+1,list.begin()+list.size()-1), base); + } + std::vector points_ranges; + eckit::Tokenizer{","}(list,points_ranges); + auto tokenize_range = eckit::Tokenizer{"-"}; + auto to_int = [](const std::string& s, const eckit::CodeLocation& here) -> long { + try { + return std::stol(s); + } + catch( std::exception& e) { + throw_Exception("Could not convert '" + s + "' to integer",here); + } + }; + for (const auto& point_range: points_ranges) { + std::vector point_range_split; + tokenize_range(point_range,point_range_split); + if (point_range_split.size() == 1) { + points.emplace_back(to_int(point_range,Here())-base); + } + else { + for (long p = to_int(point_range_split[0],Here()); p <= to_int(point_range_split[1],Here()); ++p) { + points.emplace_back(p-base); + } + } + } + return points; +} +} + +//------------------------------------------------------------------------------ + +GridPointsJSONWriter::GridPointsJSONWriter(Grid grid, const eckit::Parametrisation& args) : grid_{grid} { + args.get("json.precision",precision_=-1); + args.get("verbose",verbose_=0); + if (not args.get("partitions",nb_partitions_=0)) { + args.get("partitioner.partitions",nb_partitions_); + } + if (not args.get("partition",partition_=-1)) { + args.get("partition",partition_); + } + args.get("json.pretty", pretty_=false); + args.get("field",field_="lonlat"); + args.get("field_base",field_base_=0); + std::string points_list; + if (args.get("index",points_list)) { + args.get("index_base",points_base_ = 0); + points_ = points_from_list(points_list,points_base_); + } + + if( nb_partitions_ > 0 ) { + auto partitioner_config = grid_.partitioner(); // recommended by the grid itself + std::string partitioner; + if (args.get("partitioner.type",partitioner)) { + partitioner_config.set("type", partitioner); + } + partitioner_config.set("partitions", nb_partitions_); + distribution_ = grid::Distribution{grid_,partitioner_config}; + } +} + +//------------------------------------------------------------------------------ + +void GridPointsJSONWriter::write(std::ostream& out, eckit::Channel& info) const { + write(out, &info); +} + +//------------------------------------------------------------------------------ + +void GridPointsJSONWriter::write(std::ostream& out, std::ostream* info) const { + int points_newline = pretty_ ? true : false; + int points_indent = pretty_ ? 2 : 0; + int partition_indent = 0; + size_t chunk_size = 1000000; + + auto saved_fmtflags = out.flags(); + if (precision_ >= 0) { + out << std::fixed << std::setprecision(precision_); + } + + auto writePoint = [&out](const Point2& p) { + out << "[" << p[0] << "," << p[1] << "]"; + }; + + if(nb_partitions_ == 0 || field_ == "partition") { + out << "["; + if (pretty_) { + out << '\n'; + } + std::string join = points_newline ? ",\n" : ","; + std::string indent(points_indent,' '); + + size_t n{0}; + if (points_.size()) { + if (field_ == "lonlat") { + auto lonlat = grid_.lonlat(); + for (auto p: points_) { + if (p < grid_.size()) { + if (n!=0) { + out << join; + } + out << indent; + auto it = lonlat.begin()+(p); + writePoint(*it); + ++n; + } + } + } + else if (field_ == "index") { + for (auto p: points_) { + if (p < grid_.size()) { + if (n!=0) { + out << join; + } + out << indent; + out << p + field_base_; + ++n; + } + } + } + else if (field_ == "partition") { + for (auto p: points_) { + if (p < grid_.size()) { + if (n!=0) { + out << join; + } + out << indent; + out << distribution_.partition(p) + field_base_; + ++n; + } + } + } + else { + ATLAS_THROW_EXCEPTION("Cannot output field \""< chunk_size) { + *info << std::fixed << std::setprecision(1) << double(begin)/double(grid_.size()) * 100 << "% completed" << std::endl; + } + }; + auto write_progress_end = [&]() { + if (info && verbose_ && grid_.size() > chunk_size) { + *info << "100% completed" << std::endl; + } + }; + while(end != grid_.size()) { + begin = end; + end = std::min(grid_.size(),begin+chunk_size); + write_progress(); + if (field_ == "lonlat") { + auto it = grid_.lonlat().begin()+begin; + for (size_t n=begin; n part(std::min(chunk_size,grid_.size())); + auto write_chunk = [&](int partition, size_t begin, size_t end, size_t& n) { + size_t size{end - begin}; + if (partition != -1) { + distribution_.partition(begin, end, part); + auto part_begin = part.begin(); + auto part_end = part.begin()+size; + if (std::find(part_begin,part_end,partition) == part_end) { + return; // no points in this chunk to write + } + } + else { + part.assign(part.size(),-1); + } + std::string join = points_newline ? ",\n" : ","; + std::string indent(points_indent,' '); + if( field_ == "lonlat" ) { + auto it = grid_.lonlat().begin() + begin; + for (size_t j=0; j chunk_size) { + *info << std::fixed << std::setprecision(1) << double(begin)/double(grid_.size()) * 100 << "% completed" << std::endl; + } + end = std::min(grid_.size(),begin+chunk_size); + write_chunk(partition,begin,end,n); + } + if (verbose_) { + if (info && grid_.size() > chunk_size) { + *info << "100% completed" << std::endl; + } + } + if (pretty_) { + out << '\n'; + } + out << std::string(partition_indent,' ') << "]" << std::flush; + return n; + }; + + if (partition_ >= 0) { + auto n = write_partition(partition_); + out << std::endl; + if (info) { + *info << "Partition " << partition_ << " contains " << n << " points." << std::endl; + } + } + else { + if (pretty_) { + partition_indent += 2; + points_indent += 2; + } + out << "[\n"; + std::vector points_per_partition(nb_partitions_); + for( int p = 0; p < nb_partitions_; ++p ) { + if (p != 0) { + out << ",\n"; + } + points_per_partition[p] = write_partition(p); + } + out << "\n]"; + out << std::endl; + if (info) { + *info << "Points per partition:" << std::endl; + for (size_t p=0; p + +#include "eckit/config/Parametrisation.h" + +#include "atlas/grid.h" + +namespace atlas { +namespace util { + +//------------------------------------------------------------------------------ + +class GridPointsJSONWriter { +public: + GridPointsJSONWriter(Grid grid, const eckit::Parametrisation& args); + + void write(std::ostream& out, eckit::Channel& info) const; + + void write(std::ostream& out, std::ostream* info = nullptr) const; + +private: + Grid grid_; + grid::Distribution distribution_; + int precision_; + int verbose_; + int nb_partitions_; + int partition_; + bool pretty_; + std::string field_; + std::vector points_; + long points_base_; + long field_base_; +}; + +//------------------------------------------------------------------------------ + +} // namespace util +} // namespace atlas \ No newline at end of file From 752d82f81bb3207a371a6251bb1220a43501e965 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Oct 2023 11:37:08 +0200 Subject: [PATCH 21/24] Better error message when MatchingPartitioner is not initialised with Mesh/FunctionSpace to match --- .../detail/partitioner/MatchingFunctionSpacePartitioner.cc | 4 +++- src/atlas/grid/detail/partitioner/MatchingMeshPartitioner.cc | 5 ++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/atlas/grid/detail/partitioner/MatchingFunctionSpacePartitioner.cc b/src/atlas/grid/detail/partitioner/MatchingFunctionSpacePartitioner.cc index 44546a090..601a935d7 100644 --- a/src/atlas/grid/detail/partitioner/MatchingFunctionSpacePartitioner.cc +++ b/src/atlas/grid/detail/partitioner/MatchingFunctionSpacePartitioner.cc @@ -10,6 +10,8 @@ #include "atlas/grid/detail/partitioner/MatchingFunctionSpacePartitioner.h" +#include + #include "atlas/functionspace/FunctionSpace.h" #include "atlas/runtime/Exception.h" @@ -19,7 +21,7 @@ namespace detail { namespace partitioner { MatchingFunctionSpacePartitioner::MatchingFunctionSpacePartitioner(): Partitioner() { - ATLAS_NOTIMPLEMENTED; + ATLAS_THROW_EXCEPTION("Error: A MatchingFunctionPartitioner needs to be initialised with a FunctionSpace"); } // MatchingFunctionSpacePartitioner::MatchingFunctionSpacePartitioner(const idx_t nb_partitions): diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitioner.cc b/src/atlas/grid/detail/partitioner/MatchingMeshPartitioner.cc index 3c256c65e..6ae913afc 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitioner.cc +++ b/src/atlas/grid/detail/partitioner/MatchingMeshPartitioner.cc @@ -9,6 +9,9 @@ */ #include "atlas/grid/detail/partitioner/MatchingMeshPartitioner.h" + +#include + #include "atlas/runtime/Exception.h" namespace atlas { @@ -17,7 +20,7 @@ namespace detail { namespace partitioner { MatchingMeshPartitioner::MatchingMeshPartitioner(): Partitioner() { - ATLAS_NOTIMPLEMENTED; + ATLAS_THROW_EXCEPTION("Error: A MatchingMeshPartitioner needs to be initialised with a Mesh"); } // MatchingMeshPartitioner::MatchingMeshPartitioner(const idx_t nb_partitions): Partitioner(nb_partitions) { From 2c558fcf48275d4d5f1610f29b2c19502e9854d2 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Oct 2023 15:02:14 +0200 Subject: [PATCH 22/24] Add errors for healpix grid when ordering is not 'ring' --- src/atlas/grid/StructuredGrid.cc | 2 +- src/atlas/grid/StructuredGrid.h | 2 +- src/atlas/grid/detail/grid/Healpix.cc | 20 +++++++++++++------- src/atlas/grid/detail/grid/Healpix.h | 2 +- src/tests/mesh/test_healpixmeshgen.cc | 20 ++++++++++++++++++++ 5 files changed, 36 insertions(+), 10 deletions(-) diff --git a/src/atlas/grid/StructuredGrid.cc b/src/atlas/grid/StructuredGrid.cc index 86c0c389d..7ade36468 100644 --- a/src/atlas/grid/StructuredGrid.cc +++ b/src/atlas/grid/StructuredGrid.cc @@ -76,6 +76,6 @@ inline const HealpixGrid::grid_t* healpix_grid(const Grid::Implementation* grid) HealpixGrid::HealpixGrid(const Grid& grid): StructuredGrid(grid), grid_(healpix_grid(get())) {} -HealpixGrid::HealpixGrid(int N): HealpixGrid(Grid(new HealpixGrid::grid_t(N))) {} +HealpixGrid::HealpixGrid(int N, const std::string& ordering): HealpixGrid(Grid(new HealpixGrid::grid_t(N, ordering))) {} } // namespace atlas diff --git a/src/atlas/grid/StructuredGrid.h b/src/atlas/grid/StructuredGrid.h index a78feeabb..16d3893c7 100644 --- a/src/atlas/grid/StructuredGrid.h +++ b/src/atlas/grid/StructuredGrid.h @@ -272,7 +272,7 @@ class HealpixGrid : public StructuredGrid { public: HealpixGrid(const Grid&); - HealpixGrid(int N); + HealpixGrid(int N, const std::string& ordering = "ring"); operator bool() const { return valid(); } diff --git a/src/atlas/grid/detail/grid/Healpix.cc b/src/atlas/grid/detail/grid/Healpix.cc index 9a0d4c6f3..1150b54db 100644 --- a/src/atlas/grid/detail/grid/Healpix.cc +++ b/src/atlas/grid/detail/grid/Healpix.cc @@ -14,6 +14,7 @@ #include #include #include +#include #include "atlas/grid/Spacing.h" #include "atlas/grid/StructuredGrid.h" @@ -42,16 +43,17 @@ static class HealpixGridBuilder : GridBuilder { int id; std::vector matches; if (match(name, matches, id)) { - int N = std::stoi(matches[0]); - return new detail::grid::Healpix(N); + int N = std::stoi(matches[0]); + std::string ordering = config.getString("ordering", "ring"); + return new detail::grid::Healpix(N, ordering); } return nullptr; } const Grid::Implementation* create(const Grid::Config& config) const override { - long N; - config.get("N", N); - return new detail::grid::Healpix(N); + long N = config.getLong("N", 0); + std::string ordering = config.getString("ordering", "ring"); + return new detail::grid::Healpix(N, ordering); } } healpix_builder_; @@ -110,8 +112,12 @@ Healpix::YSpace healpix_yspace(long N) { return new spacing::CustomSpacing(y.size(), y.data()); } -Healpix::Healpix(long N): - Structured("H" + std::to_string(N), healpix_xspace(N), healpix_yspace(N), Projection(), GlobalDomain()) {} +Healpix::Healpix(long N, const std::string& ordering): + Structured("H" + std::to_string(N), healpix_xspace(N), healpix_yspace(N), Projection(), GlobalDomain()) { + if (ordering != "ring") { + ATLAS_THROW_EXCEPTION("atlas Healpix Grid is only supported with ring ordering"); + } + } Healpix::Config Healpix::meshgenerator() const { diff --git a/src/atlas/grid/detail/grid/Healpix.h b/src/atlas/grid/detail/grid/Healpix.h index fbe2de34a..3d5f4dd92 100644 --- a/src/atlas/grid/detail/grid/Healpix.h +++ b/src/atlas/grid/detail/grid/Healpix.h @@ -20,7 +20,7 @@ namespace grid { class Healpix : public Structured { public: using Structured::Structured; - Healpix(long N); + Healpix(long N, const std::string& ordering = "ring"); Config meshgenerator() const override; Config partitioner() const override; static std::string static_type() { return "healpix"; } diff --git a/src/tests/mesh/test_healpixmeshgen.cc b/src/tests/mesh/test_healpixmeshgen.cc index c83be6aaa..e49bdba8f 100644 --- a/src/tests/mesh/test_healpixmeshgen.cc +++ b/src/tests/mesh/test_healpixmeshgen.cc @@ -128,10 +128,30 @@ CASE("test_create_healpix_mesh") { //----------------------------------------------------------------------------- +CASE("construction by config") { + EXPECT(Grid(util::Config("type", "healpix")("N", 3)) == Grid("H3")); +} + CASE("construction by integer") { EXPECT(HealpixGrid(3) == Grid("H3")); } +CASE("construction by integer and ordering") { + EXPECT(HealpixGrid(3, "ring") == Grid("H3")); +} + +CASE("failing construction by integer and unsupported ordering") { + EXPECT_THROWS(HealpixGrid grid(3, "nested")); +} + +CASE("failing construction by config with unsupported ordering") { + EXPECT_THROWS(Grid grid(util::Config("name", "H3")("ordering", "nested"))); +} + +CASE("failing construction by config with unsupported ordering") { + EXPECT_THROWS(Grid grid(util::Config("type", "healpix")("N", 3)("ordering", "nested"))); +} + //----------------------------------------------------------------------------- CASE("matching mesh partitioner") { From eb61496d044fe13b01c1b3ea8533f126e6182402 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Oct 2023 15:05:41 +0200 Subject: [PATCH 23/24] Version 0.35.1 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 7b52f5e51..731b95d7f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.35.0 +0.35.1 From 7e6c7d24dd89f4b2d7b49dbc44f457adcb0baa74 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 24 Oct 2023 15:38:13 +0200 Subject: [PATCH 24/24] Update changelog --- CHANGELOG.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7630786a6..15d1a9c05 100644 --- 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 + +### 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