From e734e7a9a36417d0132e4a030069fd47bab6f0e0 Mon Sep 17 00:00:00 2001 From: Oliver Lomax Date: Fri, 15 Dec 2023 13:49:13 +0000 Subject: [PATCH 01/23] Add SphericalVector interpolation method using parallel transport (#163) * Implemented forward interpolation. TODO: adjoint and more tests. * Fixed race condition errors. * Replaced rotation matrices with complex interpolation weights. * Fixed typos before discussion. * Renamed class to SphericalVector. Added Eigen3 matrices. * Update SphericalVector.cc make_view calls erroneously in loop body. * Tidied up SphericalVector class. * Minor cosmetic changes. * Replaced optional compilation with #if ATLAS_HAVE_EIGEN in source file and header. * Removed redundant macros. * Removed static factory linking. * Fused horizontal and vertical component matrix-multiplications. TODO: Write test for 3-vector field (note: this works when hacking current 2-vector test). * Tidied fused loop. * Uncovered and fixed differences in eckit and Eigen3 CRS format. Also added tests for 3-vector and 2 vector fields. * Added multiple levels to 3d fields. * Add SphericalVector to MethodFactory * Added more consistent types to iteration indices. * Further index consistency added. * Removed superfluous templates. * Tided up macros. * Disable OpenMP for older intel-classic compiler (< intel/2022.2) * Enable test with CONDITION statement * Revert whitespace changes * Make greatCircleCourse private before moving to eckit * Fix header includes * Addressed reviewer comments. --- src/atlas/CMakeLists.txt | 2 + src/atlas/functionspace/NodeColumns.cc | 4 + src/atlas/interpolation/method/Method.h | 4 +- .../interpolation/method/MethodFactory.cc | 2 + .../method/sphericalvector/SphericalVector.cc | 317 ++++++++++++++++++ .../method/sphericalvector/SphericalVector.h | 120 +++++++ src/atlas/util/Geometry.cc | 44 ++- src/atlas/util/Geometry.h | 23 ++ src/tests/interpolation/CMakeLists.txt | 10 +- .../test_interpolation_cubedsphere.cc | 6 +- .../test_interpolation_spherical_vector.cc | 241 +++++++++++++ src/tests/util/CMakeLists.txt | 6 + src/tests/util/test_unitsphere.cc | 43 +++ 13 files changed, 815 insertions(+), 7 deletions(-) create mode 100644 src/atlas/interpolation/method/sphericalvector/SphericalVector.cc create mode 100644 src/atlas/interpolation/method/sphericalvector/SphericalVector.h create mode 100644 src/tests/interpolation/test_interpolation_spherical_vector.cc create mode 100644 src/tests/util/test_unitsphere.cc diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 02d2ec38d..3636d74d6 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -634,6 +634,8 @@ interpolation/method/knn/KNearestNeighboursBase.cc interpolation/method/knn/KNearestNeighboursBase.h interpolation/method/knn/NearestNeighbour.cc interpolation/method/knn/NearestNeighbour.h +interpolation/method/sphericalvector/SphericalVector.h +interpolation/method/sphericalvector/SphericalVector.cc interpolation/method/structured/Cubic2D.cc interpolation/method/structured/Cubic2D.h interpolation/method/structured/Cubic3D.cc diff --git a/src/atlas/functionspace/NodeColumns.cc b/src/atlas/functionspace/NodeColumns.cc index 86259313b..76fe397e3 100644 --- a/src/atlas/functionspace/NodeColumns.cc +++ b/src/atlas/functionspace/NodeColumns.cc @@ -283,6 +283,10 @@ void NodeColumns::set_field_metadata(const eckit::Configuration& config, Field& idx_t variables(0); config.get("variables", variables); field.set_variables(variables); + + if (config.has("type")) { + field.metadata().set("type", config.getString("type")); + } } array::DataType NodeColumns::config_datatype(const eckit::Configuration& config) const { diff --git a/src/atlas/interpolation/method/Method.h b/src/atlas/interpolation/method/Method.h index 43bd1f4d8..5e5abbfde 100644 --- a/src/atlas/interpolation/method/Method.h +++ b/src/atlas/interpolation/method/Method.h @@ -127,6 +127,8 @@ class Method : public util::Object { virtual void do_setup(const FunctionSpace& source, const Field& target); virtual void do_setup(const FunctionSpace& source, const FieldSet& target); + void check_compatibility(const Field& src, const Field& tgt, const Matrix& W) const; + private: template void interpolate_field(const Field& src, Field& tgt, const Matrix&) const; @@ -152,8 +154,6 @@ class Method : public util::Object { template void adjoint_interpolate_field_rank3(Field& src, const Field& tgt, const Matrix&) const; - void check_compatibility(const Field& src, const Field& tgt, const Matrix& W) const; - private: const Matrix* matrix_ = nullptr; std::shared_ptr matrix_shared_; diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index 7c9f6952c..4baa44c37 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -16,6 +16,7 @@ #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" #include "knn/NearestNeighbour.h" +#include "sphericalvector/SphericalVector.h" #include "structured/Cubic2D.h" #include "structured/Cubic3D.h" #include "structured/Linear2D.h" @@ -47,6 +48,7 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); + MethodBuilder(); } } link; } diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc new file mode 100644 index 000000000..483ce8634 --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -0,0 +1,317 @@ +/* + * (C) Crown Copyright 2023 Met Office + * + * 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. + */ + +#include "atlas/library/defines.h" +#include "atlas/interpolation/method/sphericalvector/SphericalVector.h" + +#include +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/array/ArrayView.h" +#include "atlas/array/helpers/ArrayForEach.h" +#include "atlas/array/Range.h" +#include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" +#include "atlas/interpolation/Cache.h" +#include "atlas/interpolation/Interpolation.h" +#include "atlas/interpolation/method/MethodFactory.h" +#include "atlas/parallel/omp/omp.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Trace.h" +#include "atlas/util/Constants.h" +#include "atlas/util/Geometry.h" + +namespace atlas { +namespace interpolation { +namespace method { + +namespace { +MethodBuilder __builder("spherical-vector"); +} + +#if ATLAS_HAVE_EIGEN + +// A bug exists in intel versions < intel/2022.2 with OpenMP +// Intel OneAPI version 2022.2 corresponds to Intel classic (icpc) version 2021.6 +#if defined(__INTEL_COMPILER) && defined(__INTEL_COMPILER_UPDATE) +#if (__INTEL_COMPILER <= 2021) && (__INTEL_COMPILER_UPDATE < 6) +#warning Disabling OpenMP to prevent internal compiler error for intel-classic version < 2021.6 (intel-oneapi/2022.2) +#undef atlas_omp_parallel_for +#define atlas_omp_parallel_for for +#endif +#endif + +using Complex = SphericalVector::Complex; + +template +using SparseMatrix = SphericalVector::SparseMatrix; +using RealMatrixMap = Eigen::Map>; +using ComplexTriplets = std::vector>; +using RealTriplets = std::vector>; + +using EckitMatrix = eckit::linalg::SparseMatrix; + +namespace { + +RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { + return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(), + baseMatrix.nonZeros(), baseMatrix.outer(), + baseMatrix.inner(), baseMatrix.data()); +} + +template +auto getInnerIt(const Matrix& matrix, typename Matrix::Index k) { + return typename Matrix::InnerIterator(matrix, k); +} + +template +void sparseMatrixForEach(const Functor& functor, const Matrix& matrix) { + using Index = typename Matrix::Index; + atlas_omp_parallel_for (auto k = Index{}; k < matrix.outerSize(); ++k) { + for (auto it = getInnerIt(matrix, k); it; ++it) { + functor(it.row(), it.col(), it.value()); + } + } +} + +template +void sparseMatrixForEach(const Functor& functor, const Matrix1& matrix1, + const Matrix2& matrix2) { + using Index = typename Matrix1::Index; + atlas_omp_parallel_for (auto k = Index{}; k < matrix1.outerSize(); ++k) { + for (auto [it1, it2] = + std::make_pair(getInnerIt(matrix1, k), getInnerIt(matrix2, k)); + it1; ++it1, ++it2) { + functor(it1.row(), it1.col(), it1.value(), it2.value()); + } + } +} + +template +void matrixMultiply(const SourceView& sourceView, TargetView& targetView, + const Functor& multiplyFunctor, + const Matrices&... matrices) { + + const auto multiplyColumn = [&](auto i, auto j, const auto&... weights) { + constexpr auto Rank = std::decay_t::rank(); + if constexpr (Rank == 2) { + const auto sourceSlice = sourceView.slice(j, array::Range::all()); + auto targetSlice = targetView.slice(i, array::Range::all()); + multiplyFunctor(sourceSlice, targetSlice, weights...); + } else if constexpr(Rank == 3) { + const auto sourceSlice = + sourceView.slice(j, array::Range::all(), array::Range::all()); + auto targetSlice = + targetView.slice(i, array::Range::all(), array::Range::all()); + array::helpers::ArrayForEach<0>::apply( + std::tie(sourceSlice, targetSlice), + [&](auto&& sourceVars, auto&& targetVars) { + multiplyFunctor(sourceVars, targetVars, weights...); + }); + } else { + ATLAS_NOTIMPLEMENTED; + } + }; + + sparseMatrixForEach(multiplyColumn, matrices...); +} + +} // namespace + +void SphericalVector::do_setup(const Grid& source, const Grid& target, + const Cache&) { + ATLAS_NOTIMPLEMENTED; +} + +void SphericalVector::do_setup(const FunctionSpace& source, + const FunctionSpace& target) { + ATLAS_TRACE("interpolation::method::SphericalVector::do_setup"); + source_ = source; + target_ = target; + + if (target_.size() == 0) { + return; + } + + setMatrix(Interpolation(interpolationScheme_, source_, target_)); + + // Get matrix data. + const auto nRows = matrix().rows(); + const auto nCols = matrix().cols(); + const auto nNonZeros = matrix().nonZeros(); + const auto baseWeights = makeMatrixMap(matrix()); + + // Note: need to store copy of weights as Eigen3 sorts compressed rows by j + // whereas eckit does not. + complexWeights_ = std::make_shared(nRows, nCols); + realWeights_ = std::make_shared(nRows, nCols); + auto complexTriplets = ComplexTriplets(nNonZeros); + auto realTriplets = RealTriplets(nNonZeros); + + const auto sourceLonLats = array::make_view(source_.lonlat()); + const auto targetLonLats = array::make_view(target_.lonlat()); + + geometry::UnitSphere unitSphere; + + const auto setWeights = [&](auto i, auto j, const auto& baseWeight) { + const auto sourceLonLat = + PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); + const auto targetLonLat = + PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); + + const auto alpha = unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); + + const auto deltaAlpha = + (alpha.first - alpha.second) * util::Constants::degreesToRadians(); + + const auto idx = std::distance(baseWeights.valuePtr(), &baseWeight); + + complexTriplets[idx] = {int(i), int(j), std::polar(baseWeight, deltaAlpha)}; + realTriplets[idx] = {int(i), int(j), baseWeight}; + }; + + sparseMatrixForEach(setWeights, baseWeights); + complexWeights_->setFromTriplets(complexTriplets.begin(), + complexTriplets.end()); + realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end()); + + ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); + ATLAS_ASSERT(realWeights_->nonZeros() == matrix().nonZeros()); +} + +SphericalVector::SphericalVector(const Config& config) : Method(config) { + const auto& conf = dynamic_cast(config); + interpolationScheme_ = conf.getSubConfiguration("scheme"); +} + +void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } + +void SphericalVector::do_execute(const FieldSet& sourceFieldSet, + FieldSet& targetFieldSet, + Metadata& metadata) const { + ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); + + const auto nFields = sourceFieldSet.size(); + ATLAS_ASSERT(nFields == targetFieldSet.size()); + + for (auto i = 0; i < sourceFieldSet.size(); ++i) { + do_execute(sourceFieldSet[i], targetFieldSet[i], metadata); + } +} + +void SphericalVector::do_execute(const Field& sourceField, Field& targetField, + Metadata&) const { + ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); + + const auto fieldType = sourceField.metadata().getString("type", ""); + if (fieldType != "vector") { + + auto metadata = Metadata(); + Method::do_execute(sourceField, targetField, metadata); + + return; + } + + if (target_.size() == 0) { + return; + } + + ATLAS_ASSERT_MSG(sourceField.variables() == 2 || sourceField.variables() == 3, + "Vector field can only have 2 or 3 components."); + + Method::check_compatibility(sourceField, targetField, matrix()); + + haloExchange(sourceField); + + if (sourceField.datatype().kind() == array::DataType::KIND_REAL64) { + interpolate_vector_field(sourceField, targetField); + } else if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { + interpolate_vector_field(sourceField, targetField); + } else { + ATLAS_NOTIMPLEMENTED; + } + + targetField.set_dirty(); +} + +void SphericalVector::do_execute_adjoint(FieldSet& sourceFieldSet, + const FieldSet& targetFieldSet, + Metadata& metadata) const { + ATLAS_NOTIMPLEMENTED; +} + +void SphericalVector::do_execute_adjoint(Field& sourceField, + const Field& targetField, + Metadata& metadata) const { + ATLAS_NOTIMPLEMENTED; +} + +template +void SphericalVector::interpolate_vector_field(const Field& sourceField, + Field& targetField) const { + if (sourceField.rank() == 2) { + interpolate_vector_field(sourceField, targetField); + return; + } + + if (sourceField.rank() == 3) { + interpolate_vector_field(sourceField, targetField); + return; + } + + ATLAS_NOTIMPLEMENTED; +} + +template +void SphericalVector::interpolate_vector_field(const Field& sourceField, + Field& targetField) const { + + const auto sourceView = array::make_view(sourceField); + auto targetView = array::make_view(targetField); + targetView.assign(0.); + + const auto horizontalComponent = [](const auto& sourceVars, auto& targetVars, + const auto& complexWeight) { + const auto sourceVector = Complex(sourceVars(0), sourceVars(1)); + const auto targetVector = complexWeight * sourceVector; + targetVars(0) += targetVector.real(); + targetVars(1) += targetVector.imag(); + }; + + if (sourceField.variables() == 2) { + matrixMultiply(sourceView, targetView, horizontalComponent, + *complexWeights_); + return; + } + + if (sourceField.variables() == 3) { + + const auto horizontalAndVerticalComponent = [&]( + const auto& sourceVars, auto& targetVars, const auto& complexWeight, + const auto& realWeight) { + horizontalComponent(sourceVars, targetVars, complexWeight); + targetVars(2) += realWeight * sourceVars(2); + }; + + matrixMultiply(sourceView, targetView, horizontalAndVerticalComponent, + *complexWeights_, *realWeights_); + return; + } + + ATLAS_NOTIMPLEMENTED; +} + +#endif + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h new file mode 100644 index 000000000..91ce2a1a2 --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -0,0 +1,120 @@ +/* + * (C) Crown Copyright 2023 Met Office + * + * 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. + */ + +#include "atlas/library/defines.h" + +#pragma once + +#include +#include + +#if ATLAS_HAVE_EIGEN +#include +#endif + +#include "atlas/functionspace/FunctionSpace.h" +#include "atlas/interpolation/method/Method.h" +#include "atlas/linalg/sparse.h" + + +namespace atlas { +namespace interpolation { +namespace method { + +#if ATLAS_HAVE_EIGEN +class SphericalVector : public Method { + public: + using Complex = std::complex; + + template + using SparseMatrix = Eigen::SparseMatrix; + using ComplexMatrix = SparseMatrix; + using RealMatrix = SparseMatrix; + + + /// @brief Interpolation post-processor for vector field interpolation + /// + /// @details Takes a base interpolation config keyed to "scheme" and creates + /// A set of complex intepolation weights which rotate source vector + /// elements into the target elements' individual frame of reference. + /// Method works by creating a great-circle arc between the source + /// and target locations; the amount of rotation is determined by the + /// difference the in the great-cricle course (relative to north) at + /// the source and target location. + /// Both source and target fields require the "type" metadata to be + /// set to "vector" for this method to be invoked. Otherwise, the + /// base scalar interpolation method is invoked. + /// Note: This method only works with matrix-based interpolation + /// schemes. + /// + SphericalVector(const Config& config); + ~SphericalVector() override {} + + void print(std::ostream&) const override; + const FunctionSpace& source() const override { return source_; } + const FunctionSpace& target() const override { return target_; } + + void do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, + Metadata& metadata) const override; + void do_execute(const Field& sourceField, Field& targetField, + Metadata& metadata) const override; + + void do_execute_adjoint(FieldSet& sourceFieldSet, + const FieldSet& targetFieldSet, + Metadata& metadata) const override; + void do_execute_adjoint(Field& sourceField, const Field& targetField, + Metadata& metadata) const override; + + private: + template + void interpolate_vector_field(const Field& source, Field& target) const; + + template + void interpolate_vector_field(const Field& source, Field& target) const; + + using Method::do_setup; + void do_setup(const FunctionSpace& source, + const FunctionSpace& target) override; + void do_setup(const Grid& source, const Grid& target, const Cache&) override; + + eckit::LocalConfiguration interpolationScheme_; + + FunctionSpace source_; + FunctionSpace target_; + + std::shared_ptr complexWeights_; + std::shared_ptr realWeights_; + +}; +#else + class SphericalVector : public Method { + public: + SphericalVector(const Config& config) : Method(config) { + ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen"); + } + + ~SphericalVector() override {} + + void print(std::ostream&) const override {} + const FunctionSpace& source() const override {ATLAS_NOTIMPLEMENTED;} + const FunctionSpace& target() const override {ATLAS_NOTIMPLEMENTED;} + + void do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, + Metadata& metadata) const override {} + void do_execute(const Field& sourceField, Field& targetField, + Metadata& metadata) const override {} + private: + void do_setup(const FunctionSpace& source, + const FunctionSpace& target) override {} + void do_setup(const Grid& source, const Grid& target, const Cache&) override {} + }; +#endif + + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/util/Geometry.cc b/src/atlas/util/Geometry.cc index 7d8e07797..074d72d34 100644 --- a/src/atlas/util/Geometry.cc +++ b/src/atlas/util/Geometry.cc @@ -8,12 +8,16 @@ * nor does it submit to any jurisdiction. */ +#include "atlas/util/Geometry.h" + +#include + #include "eckit/geometry/Point2.h" #include "eckit/geometry/Point3.h" #include "atlas/library/config.h" #include "atlas/runtime/Exception.h" -#include "atlas/util/Geometry.h" +#include "atlas/util/Constants.h" namespace atlas { @@ -30,6 +34,44 @@ void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const { Sphere::convertCartesianToSpherical(radius_, xyz, lonlat); } + +/// @brief Calculate great-cricle course between points +/// +/// @details Calculates the direction (clockwise from north) of a great-circle +/// arc between lonLat1 and lonLat2. Returns the direction of the arc +/// at lonLat1 (first) and lonLat2 (second). Angle is normalised to the +/// range of atan2 (usually (-180, 180]). All input and output values +/// are in units of degrees. +/// @ref https://en.wikipedia.org/wiki/Great-circle_navigation +/// +std::pair greatCircleCourse(const Point2& lonLat1, + const Point2& lonLat2) { + + const auto lambda1 = lonLat1[0] * util::Constants::degreesToRadians(); + const auto lambda2 = lonLat2[0] * util::Constants::degreesToRadians(); + const auto phi1 = lonLat1[1] * util::Constants::degreesToRadians(); + const auto phi2 = lonLat2[1] * util::Constants::degreesToRadians(); + + const auto sinLambda12 = std::sin(lambda2 - lambda1); + const auto cosLambda12 = std::cos(lambda2 - lambda1); + const auto sinPhi1 = std::sin(phi1); + const auto sinPhi2 = std::sin(phi2); + const auto cosPhi1 = std::cos(phi1); + const auto cosPhi2 = std::cos(phi2); + + const auto alpha1 = + std::atan2(cosPhi2 * sinLambda12, + cosPhi1 * sinPhi2 - sinPhi1 * cosPhi2 * cosLambda12); + + const auto alpha2 = + std::atan2(cosPhi1 * sinLambda12, + -cosPhi2 * sinPhi1 + sinPhi2 * cosPhi1 * cosLambda12); + + return std::make_pair(alpha1 * util::Constants::radiansToDegrees(), + alpha2 * util::Constants::radiansToDegrees()); +}; + + } // namespace detail } // namespace geometry diff --git a/src/atlas/util/Geometry.h b/src/atlas/util/Geometry.h index 3880645e8..dd8173c3a 100644 --- a/src/atlas/util/Geometry.h +++ b/src/atlas/util/Geometry.h @@ -13,6 +13,7 @@ #include #include #include +#include #include "atlas/runtime/Exception.h" #include "atlas/util/Earth.h" @@ -27,6 +28,20 @@ namespace geometry { namespace detail { +// TODO: move greatCircleCourse to eckit::geometry::Sphere + +/// @brief Calculate great-cricle course between points +/// +/// @details Calculates the direction (clockwise from north) of a great-circle +/// arc between lonLat1 and lonLat2. Returns the direction of the arc +/// at lonLat1 (first) and lonLat2 (second). Angle is normalised to the +/// range of atan2 (usually (-180, 180]). All input and output values +/// are in units of degrees. +/// @ref https://en.wikipedia.org/wiki/Great-circle_navigation +std::pair greatCircleCourse(const Point2& lonLat1, const Point2& lonLat2); + +//------------------------------------------------------------------------------------------------------ + class GeometryBase : public util::Object { public: virtual ~GeometryBase() = default; @@ -36,6 +51,7 @@ class GeometryBase : public util::Object { virtual double distance(const Point3& p1, const Point3& p2) const = 0; virtual double radius() const = 0; virtual double area() const = 0; + virtual std::pair greatCircleCourse(const Point2& p1, const Point2& p2) = 0; Point3 xyz(const Point2& lonlat) const { Point3 xyz; @@ -64,6 +80,9 @@ class GeometrySphereT : public GeometryBase { double distance(const Point3& p1, const Point3& p2) const override { return SphereT::distance(p1, p2); } double radius() const override { return SphereT::radius(); } double area() const override { return SphereT::area(); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) override { + return atlas::geometry::detail::greatCircleCourse(p1, p2); + } }; class GeometrySphere : public GeometryBase { @@ -77,6 +96,9 @@ class GeometrySphere : public GeometryBase { double distance(const Point3& p1, const Point3& p2) const override { return Sphere::distance(radius_, p1, p2); } double radius() const override { return radius_; } double area() const override { return Sphere::area(radius_); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) override { + return atlas::geometry::detail::greatCircleCourse(p1, p2); + } private: double radius_; @@ -107,6 +129,7 @@ class Geometry : public util::ObjectHandle { double distance(const Point3& p1, const Point3& p2) const { return get()->distance(p1, p2); } double radius() const { return get()->radius(); } double area() const { return get()->area(); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) { return get()->greatCircleCourse(p1, p2); } protected: template diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 596729209..29d220590 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -122,4 +122,12 @@ ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -endif() \ No newline at end of file +ecbuild_add_test( TARGET atlas_test_interpolation_spherical_vector + OMP 4 + SOURCES test_interpolation_spherical_vector.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + CONDITION atlas_HAVE_EIGEN +) + +endif() diff --git a/src/tests/interpolation/test_interpolation_cubedsphere.cc b/src/tests/interpolation/test_interpolation_cubedsphere.cc index b2d446ebf..a6d930983 100644 --- a/src/tests/interpolation/test_interpolation_cubedsphere.cc +++ b/src/tests/interpolation/test_interpolation_cubedsphere.cc @@ -60,7 +60,7 @@ void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { // Return (u, v) field with vortex_rollup as the streamfunction. // This has no physical significance, but it makes a nice swirly field. -std::pair vortexField(double lon, double lat) { +std::pair vortexHorizontal(double lon, double lat) { // set hLon and hLat step size. const double hLon = 0.0001; @@ -222,7 +222,7 @@ CASE("cubedsphere_wind_interpolation") { const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT)); // Set (u, v) wind - std::tie(u(idx), v(idx)) = vortexField(ll.lon(), ll.lat()); + std::tie(u(idx), v(idx)) = vortexHorizontal(ll.lon(), ll.lat()); // Get wind transform jacobian. auto jac = windTransform(ll, t); @@ -282,7 +282,7 @@ CASE("cubedsphere_wind_interpolation") { std::tie(u(idx), v(idx)) = matMul(jac, vAlpha(idx), vBeta(idx)); // Get error. - const auto uvTarget = vortexField(ll.lon(), ll.lat()); + const auto uvTarget = vortexHorizontal(ll.lon(), ll.lat()); error0(idx) = Point2::distance(Point2(uvTarget.first, uvTarget.second), Point2(uOrig(idx), vOrig(idx))); error1(idx) = Point2::distance(Point2(uvTarget.first, uvTarget.second), Point2(u(idx), v(idx))); diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc new file mode 100644 index 000000000..1e68b8e91 --- /dev/null +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -0,0 +1,241 @@ +/* + * (C) Crown Copyright 2023 Met Office + * + * 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. + */ + +#include "atlas/array.h" +#include "atlas/array/helpers/ArrayForEach.h" +#include "atlas/field.h" +#include "atlas/functionspace.h" +#include "atlas/grid.h" +#include "atlas/interpolation.h" +#include "atlas/mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/util/Point.h" +#include "atlas/util/Config.h" +#include "atlas/util/Constants.h" +#include "atlas/util/function/VortexRollup.h" + +#include "tests/AtlasTestEnvironment.h" + +namespace atlas { +namespace test { + +using namespace atlas::util; +using namespace atlas::array::helpers; + +constexpr auto Rank2dField = 2; +constexpr auto Rank3dField = 3; + +// Return (u, v) field with vortex_rollup as the streamfunction. +// This has no physical significance, but it makes a nice swirly field. +std::pair vortexHorizontal(double lon, double lat) { + + // set hLon and hLat step size. + const double hLon = 0.0001; + const double hLat = 0.0001; + + // Get finite differences. + + // Set u. + const double u = (function::vortex_rollup(lon, lat + 0.5 * hLat, 0.1) - + function::vortex_rollup(lon, lat - 0.5 * hLat, 0.1)) / + hLat; + + const double v = -(function::vortex_rollup(lon + 0.5 * hLon, lat, 0.1) - + function::vortex_rollup(lon - 0.5 * hLon, lat, 0.1)) / + (hLon * std::cos(lat * util::Constants::degreesToRadians())); + + return std::make_pair(u, v); +} + +double vortexVertical(double lon, double lat) { + return function::vortex_rollup(lon, lat, 0.1); +} + +void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { + + const auto& functionSpace = fieldSet[0].functionspace(); + const auto& mesh = functionspace::NodeColumns(functionSpace).mesh(); + + const auto gmshConfig = util::Config("coordinates", "xyz") | + util::Config("ghost", true) | + util::Config("info", true); + const auto gmsh = output::Gmsh(fileName, gmshConfig); + gmsh.write(mesh); + gmsh.write(fieldSet, functionSpace); +} + +// Helper function to generate a NodeColumns functionspace +const auto generateNodeColums(const std::string& gridName, + const std::string& meshName) { + const auto grid = Grid(gridName); + const auto mesh = MeshGenerator(meshName).generate(grid); + return functionspace::NodeColumns(mesh); +} + +// Helper struct to key different Functionspaces to strings +struct FunctionSpaceFixtures { + static const FunctionSpace& get(const std::string& fixture) { + static std::map functionSpaces = { + {"cubedsphere_mesh", + generateNodeColums("CS-LFR-48", "cubedsphere_dual")}, + {"gaussian_mesh", generateNodeColums("O48", "structured")}}; + return functionSpaces.at(fixture); + } +}; + +// Helper struct to key different grid configs to strings +struct FieldSpecFixtures { + static const Config& get(const std::string& fixture) { + static std::map fieldSpecs = { + {"2vector", option::name("test field") | option::variables(2) | + option::type("vector")}, + {"3vector", option::name("test field") | option::variables(3) | + option::type("vector")}}; + return fieldSpecs.at(fixture); + } +}; + + +template +void testInterpolation(const Config& config) { + + const auto& sourceFunctionSpace = + FunctionSpaceFixtures::get(config.getString("source_fixture")); + const auto& targetFunctionSpace = + FunctionSpaceFixtures::get(config.getString("target_fixture")); + + auto sourceFieldSet = FieldSet{}; + auto targetFieldSet = FieldSet{}; + + const auto sourceLonLat = + array::make_view(sourceFunctionSpace.lonlat()); + const auto targetLonLat = + array::make_view(targetFunctionSpace.lonlat()); + + auto fieldSpec = + FieldSpecFixtures::get(config.getString("field_spec_fixture")); + if constexpr (Rank == 3) fieldSpec.set("levels", 2); + + auto sourceField = array::make_view( + sourceFieldSet.add(sourceFunctionSpace.createField(fieldSpec))); + + auto targetField = array::make_view( + targetFieldSet.add(targetFunctionSpace.createField(fieldSpec))); + + ArrayForEach<0>::apply(std::tie(sourceLonLat, sourceField), + [](auto&& lonLat, auto&& sourceColumn) { + + const auto setElems = [&](auto&& sourceElem) { + std::tie(sourceElem(0), sourceElem(1)) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (sourceElem.size() == 3) { + sourceElem(2) = vortexVertical(lonLat(0), lonLat(1)); + } + }; + if constexpr (Rank == 2) { setElems(sourceColumn); } + else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(sourceColumn), setElems); + } + }); + sourceFieldSet.set_dirty(false); + + const auto interp = Interpolation(config.getSubConfiguration("scheme"), + sourceFunctionSpace, targetFunctionSpace); + + interp.execute(sourceFieldSet, targetFieldSet); + targetFieldSet.haloExchange(); + + auto errorFieldSpec = fieldSpec; + errorFieldSpec.remove("variables"); + + auto errorField = array::make_view(targetFieldSet.add( + targetFunctionSpace.createField(errorFieldSpec))); + + auto maxError = 0.; + ArrayForEach<0>::apply(std::tie(targetLonLat, targetField, errorField), + [&](auto&& lonLat, auto&& targetColumn, + auto&& errorColumn) { + + const auto calcError = [&](auto&& targetElem, auto&& errorElem) { + auto trueValue = std::vector(targetElem.size()); + std::tie(trueValue[0], trueValue[1]) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (targetElem.size() == 3) { + trueValue[2] = vortexVertical(lonLat(0), lonLat(1)); + } + + auto errorSqrd = 0.; + for (auto k = 0; k < targetElem.size(); ++k) { + errorSqrd += + (targetElem(k) - trueValue[k]) * (targetElem(k) - trueValue[k]); + } + + errorElem = std::sqrt(errorSqrd); + maxError = std::max(maxError, static_cast(errorElem)); + }; + + if constexpr(Rank == 2) { calcError(targetColumn, errorColumn); } + else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), calcError); + } + }); + + EXPECT_APPROX_EQ(maxError, 0., config.getDouble("tol")); + + gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); + gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); +} + +CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { + + const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear"); + const auto interpScheme = + util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); + const auto cubedSphereConf = Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("file_id", "spherical_vector_cs2") + .set("scheme", interpScheme) + .set("tol", 0.00018); + + testInterpolation((cubedSphereConf)); +} + +CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { + + const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear"); + const auto interpScheme = + util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); + const auto cubedSphereConf = Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "3vector") + .set("file_id", "spherical_vector_cs3") + .set("scheme", interpScheme) + .set("tol", 0.00096); + + testInterpolation((cubedSphereConf)); +} + +CASE("finite element vector interpolation (2d-field, 2-vector)") { + + const auto baseInterpScheme = util::Config("type", "finite-element"); + const auto interpScheme = + util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); + const auto cubedSphereConf = Config("source_fixture", "gaussian_mesh") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("file_id", "spherical_vector_cs") + .set("scheme", interpScheme) + .set("tol", 0.00015); + + testInterpolation((cubedSphereConf)); +} +} +} + +int main(int argc, char** argv) { return atlas::test::run(argc, argv); } diff --git a/src/tests/util/CMakeLists.txt b/src/tests/util/CMakeLists.txt index 4e6e0daf0..47d8b6def 100644 --- a/src/tests/util/CMakeLists.txt +++ b/src/tests/util/CMakeLists.txt @@ -88,3 +88,9 @@ ecbuild_add_test( TARGET atlas_test_convexsphericalpolygon ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +ecbuild_add_test( TARGET atlas_test_unitsphere + SOURCES test_unitsphere.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + diff --git a/src/tests/util/test_unitsphere.cc b/src/tests/util/test_unitsphere.cc new file mode 100644 index 000000000..d390334a5 --- /dev/null +++ b/src/tests/util/test_unitsphere.cc @@ -0,0 +1,43 @@ +/* + * (C) Crown Copyright 2023 Met Office + * + * 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. + */ + +#include "atlas/util/Point.h" +#include "atlas/util/Geometry.h" + +#include "tests/AtlasTestEnvironment.h" + +namespace atlas { +namespace test { + +using namespace atlas::util; + +CASE("great-circle course") { + + geometry::UnitSphere g; + + const auto point1 = PointLonLat(-71.6, -33.0); // Valparaiso + const auto point2 = PointLonLat(121.8, 31.4); // Shanghai + const auto point3 = PointLonLat(0., 89.); + const auto point4 = PointLonLat(180., 89.); + + const auto targetCourse1 = -94.41; + const auto targetCourse2 = -78.42; + const auto targetCourse3 = 0.; + const auto targetCourse4 = 180.; + + const auto[ course1, course2 ] = g.greatCircleCourse(point1, point2); + const auto[ course3, course4 ] = g.greatCircleCourse(point3, point4); + + EXPECT_APPROX_EQ(course1, targetCourse1, 0.01); + EXPECT_APPROX_EQ(course2, targetCourse2, 0.01); + EXPECT_APPROX_EQ(course3, targetCourse3, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course4), targetCourse4, 1.e-14); +} +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { return atlas::test::run(argc, argv); } From 8435757b8ccd183f6cad2afb3deda2e15e624377 Mon Sep 17 00:00:00 2001 From: Matt Shin Date: Wed, 15 Nov 2023 11:28:54 +0000 Subject: [PATCH 02/23] Fix partitioner config exception message Improve logic to eliminate config exception message. --- .../meshgenerator/detail/CubedSphereDualMeshGenerator.cc | 4 +--- src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc index 953a84ed8..9026d92e8 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc @@ -63,9 +63,7 @@ CubedSphereDualMeshGenerator::CubedSphereDualMeshGenerator(const eckit::Parametr // Get partitioner. std::string partitioner; - try { p.get("partitioner.type", partitioner); } catch( std::exception& ) {} - p.get("partitioner", partitioner); - if( partitioner.size() ) { + if (p.get("partitioner", partitioner) && partitioner.size()) { options.set("partitioner", partitioner); } } diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc index c2bd898ab..250d42ed0 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc @@ -73,9 +73,7 @@ CubedSphereMeshGenerator::CubedSphereMeshGenerator(const eckit::Parametrisation& // Get partitioner. std::string partitioner; - try { p.get("partitioner.type", partitioner); } catch( std::exception& ) {} - p.get("partitioner", partitioner); - if( partitioner.size() ) { + if (p.get("partitioner", partitioner) && partitioner.size()) { options.set("partitioner", partitioner); } } From dc63328aa0540303ac19d9e043645b10c771367e Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 18 Jan 2024 13:53:05 +0100 Subject: [PATCH 03/23] Fix for setting of nb_edges via BuildEdges action --- src/atlas/mesh/actions/BuildEdges.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/atlas/mesh/actions/BuildEdges.cc b/src/atlas/mesh/actions/BuildEdges.cc index eb6461e91..08a568a17 100644 --- a/src/atlas/mesh/actions/BuildEdges.cc +++ b/src/atlas/mesh/actions/BuildEdges.cc @@ -380,7 +380,9 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { for (int halo = 0; halo <= mesh_halo; ++halo) { edge_start = edge_end; - edge_end += (edge_halo_offsets[halo + 1] - edge_halo_offsets[halo]); + if (halo+1 < edge_halo_offsets.size()) { + edge_end += (edge_halo_offsets[halo + 1] - edge_halo_offsets[halo]); + } if (/*sort edges based on lowest node local index = */ sort_edges) { if (sorted_edge_nodes_data.empty()) { @@ -552,6 +554,11 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { ss << "nb_edges_including_halo[" << i << "]"; mesh.metadata().set(ss.str(), nb_edges_including_halo[i]); } + for (int i = max_halo + 1; i <= mesh_halo; ++i) { + std::stringstream ss; + ss << "nb_edges_including_halo[" << i << "]"; + mesh.metadata().set(ss.str(), nb_edges_including_halo.back()); + } mesh.metadata().set("built_edges_for_halo", mesh_halo); From b8f8a4dc481077fc8fc217ff630ff59e94b6965a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 18 Jan 2024 14:02:11 +0100 Subject: [PATCH 04/23] Fix for elements that might have unassigned partition via parallel Delaunay meshgenerator --- src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc b/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc index 4ab849ba0..7004d6725 100644 --- a/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc @@ -197,6 +197,9 @@ void DelaunayMeshGenerator::generate(const Grid& grid, const grid::Distribution& element_nodes_uncertainty.insert(g_node_connectivity(jelem,2)); element_uncertainty.insert(jelem); } + else if (elem_ownership != GHOST && elem_ownership != CERTAIN) { + ATLAS_NOTIMPLEMENTED; + } } // Log::info() << "element_uncertainty" << std::endl; // for( auto& jelem: element_uncertainty ) { @@ -291,6 +294,9 @@ void DelaunayMeshGenerator::generate(const Grid& grid, const grid::Distribution& else if (e2 != UNCERTAIN) { elem_part = e2; } + if (elem_part == UNCERTAIN) { + elem_part = elem_node_partition(jelem,0); + } if (elem_part == part_) { collect_element(jelem); } From 2bb229380346972967b2931078e04643cc3a0999 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 25 Jan 2024 14:48:15 +0100 Subject: [PATCH 05/23] Github Actions: Fix macOS MPI slots --- tools/install-mpi.sh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/install-mpi.sh b/tools/install-mpi.sh index 03a5e3a35..c5c1847ae 100755 --- a/tools/install-mpi.sh +++ b/tools/install-mpi.sh @@ -62,6 +62,8 @@ case "$os" in openmpi) brew ls --versions openmpi || brew install openmpi echo "localhost slots=72" >> $(brew --prefix)/etc/openmpi-default-hostfile + echo "localhost slots=72" >> $(brew --prefix)/etc/prte-default-hostfile + # workaround for open-mpi/omp#7516 echo "setting the mca gds to hash..." echo "gds = hash" >> $(brew --prefix)/etc/pmix-mca-params.conf From 5f4cfb8cbf43693555882cfb3bc6f50adf20453b Mon Sep 17 00:00:00 2001 From: Dusan Figala Date: Mon, 15 Jan 2024 13:20:51 +0100 Subject: [PATCH 06/23] Add docs build workflow Checkout current atlas source Install doxygen Install latex Install texlive-full --- .github/docs-config.yml | 14 ++++++++++++++ .github/workflows/docs.yml | 15 +++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 .github/docs-config.yml create mode 100644 .github/workflows/docs.yml diff --git a/.github/docs-config.yml b/.github/docs-config.yml new file mode 100644 index 000000000..e08272a73 --- /dev/null +++ b/.github/docs-config.yml @@ -0,0 +1,14 @@ +# note: each step is executed in own process +build-steps: + - git clone --depth 1 https://github.com/ecmwf/atlas-docs.git $RUNNER_TEMP/atlas-docs + - sudo apt install -y -q doxygen texlive-full + - | + cd $RUNNER_TEMP/atlas-docs + make PUBLIC=1 WITH_ECKIT=1 WITH_DOXYGEN=1 ATLAS_SOURCE_DIR=$GITHUB_WORKSPACE clean html + echo "DOC_BUILD_PATH=$RUNNER_TEMP/atlas-docs/build/html" >> "$GITHUB_ENV" + +hosts: + ecmwf-sites: + space: docs + name: atlas + diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 000000000..1c1df2a90 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,15 @@ +name: docs + +on: + push: + tags: + - '**' + + workflow_dispatch: ~ + +jobs: + publish: + uses: ecmwf-actions/reusable-workflows/.github/workflows/cd-docs.yml@v2 + secrets: inherit + with: + config: .github/docs-config.yml From e6858536de7e6d1ea9de148a2e5f46515e41aa07 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 25 Jan 2024 16:15:52 +0000 Subject: [PATCH 07/23] Added arrayForEachDim method. --- src/atlas/array/helpers/ArrayForEach.h | 243 +++++++++++++------------ src/tests/array/test_array_foreach.cc | 59 +++++- 2 files changed, 187 insertions(+), 115 deletions(-) diff --git a/src/atlas/array/helpers/ArrayForEach.h b/src/atlas/array/helpers/ArrayForEach.h index 38d4d6441..f20b4c80b 100644 --- a/src/atlas/array/helpers/ArrayForEach.h +++ b/src/atlas/array/helpers/ArrayForEach.h @@ -1,5 +1,5 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * 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. @@ -7,9 +7,10 @@ #pragma once +#include #include #include -#include +#include #include "atlas/array/ArrayView.h" #include "atlas/array/Range.h" @@ -22,18 +23,19 @@ namespace atlas { namespace execution { -// As in C++17 std::execution namespace. Note: unsequenced_policy is a C++20 addition. +// As in C++17 std::execution namespace. Note: unsequenced_policy is a C++20 +// addition. class sequenced_policy {}; class unsequenced_policy {}; class parallel_unsequenced_policy {}; class parallel_policy {}; -// execution policy objects as in C++ std::execution namespace. Note: unseq is a C++20 addition. -inline constexpr sequenced_policy seq{ /*unspecified*/ }; -inline constexpr parallel_policy par{ /*unspecified*/ }; -inline constexpr parallel_unsequenced_policy par_unseq{ /*unspecified*/ }; -inline constexpr unsequenced_policy unseq{ /*unspecified*/ }; - +// execution policy objects as in C++ std::execution namespace. Note: unseq is a +// C++20 addition. +inline constexpr sequenced_policy seq{/*unspecified*/}; +inline constexpr parallel_policy par{/*unspecified*/}; +inline constexpr parallel_unsequenced_policy par_unseq{/*unspecified*/}; +inline constexpr unsequenced_policy unseq{/*unspecified*/}; // Type names for execution policy (Not in C++ standard) template @@ -64,34 +66,31 @@ constexpr std::string_view policy_name(execution_policy) { // Type check for execution policy (Not in C++ standard) template -constexpr auto is_execution_policy() { return - std::is_same_v || - std::is_same_v || - std::is_same_v || - std::is_same_v; +constexpr auto is_execution_policy() { + return std::is_same_v || + std::is_same_v || + std::is_same_v || + std::is_same_v; } template constexpr auto demote_policy() { - if constexpr(std::is_same_v) { + if constexpr (std::is_same_v) { return unseq; - } - else if constexpr (std::is_same_v) { + } else if constexpr (std::is_same_v) { return seq; - } - else { + } else { return ExecutionPolicy{}; } } template -constexpr auto is_omp_policy() { return - std::is_same_v || - std::is_same_v; +constexpr auto is_omp_policy() { + return std::is_same_v || + std::is_same_v; } - -template +template using demote_policy_t = decltype(demote_policy()); } // namespace execution @@ -102,7 +101,8 @@ namespace option { template util::Config execution_policy() { - return util::Config("execution_policy", execution::policy_name>()); + return util::Config("execution_policy", + execution::policy_name>()); } template @@ -110,7 +110,7 @@ util::Config execution_policy(T) { return execution_policy>(); } -} // namespace option +} // namespace option namespace array { namespace helpers { @@ -119,7 +119,9 @@ namespace detail { struct NoMask { template - constexpr bool operator()(Args...) const { return 0; } + constexpr bool operator()(Args...) const { + return 0; + } }; inline constexpr NoMask no_mask; @@ -131,13 +133,11 @@ constexpr auto tuplePushBack(const std::tuple& tuple, T value) { template void forEach(idx_t idxMax, const Functor& functor) { - - if constexpr(execution::is_omp_policy()) { + if constexpr (execution::is_omp_policy()) { atlas_omp_parallel_for(auto idx = idx_t{}; idx < idxMax; ++idx) { functor(idx); } - } - else { + } else { // Simple for-loop for sequenced or unsequenced execution policies. for (auto idx = idx_t{}; idx < idxMax; ++idx) { functor(idx); @@ -147,11 +147,10 @@ void forEach(idx_t idxMax, const Functor& functor) { template constexpr auto argPadding() { - if constexpr(NPad > 0) { + if constexpr (NPad > 0) { return std::tuple_cat(std::make_tuple(Range::all()), argPadding()); - } - else { + } else { return std::make_tuple(); } } @@ -159,29 +158,31 @@ constexpr auto argPadding() { template auto makeSlices(const std::tuple& slicerArgs, ArrayViewTuple&& arrayViews) { - constexpr auto nb_views = std::tuple_size_v; auto&& arrayView = std::get(arrayViews); using ArrayView = std::decay_t; - constexpr auto Dim = sizeof...(SlicerArgs); - constexpr auto Rank = ArrayView::rank(); + constexpr auto Dim = sizeof...(SlicerArgs); + constexpr auto Rank = ArrayView::rank(); - static_assert (Dim <= Rank, "Error: number of slicer arguments exceeds the rank of ArrayView."); + static_assert( + Dim <= Rank, + "Error: number of slicer arguments exceeds the rank of ArrayView."); const auto paddedArgs = std::tuple_cat(slicerArgs, argPadding()); const auto slicer = [&arrayView](const auto&... args) { return std::make_tuple(arrayView.slice(args...)); }; - if constexpr (ViewIdx == nb_views-1) { + if constexpr (ViewIdx == nb_views - 1) { return std::apply(slicer, paddedArgs); - } - else { + } else { // recurse - return std::tuple_cat(std::apply(slicer, paddedArgs), - makeSlices(slicerArgs, std::forward(arrayViews))); + return std::tuple_cat( + std::apply(slicer, paddedArgs), + makeSlices(slicerArgs, + std::forward(arrayViews))); } } @@ -192,80 +193,79 @@ template struct ArrayForEachImpl { template - static void apply(ArrayViewTuple&& arrayViews, - const Mask& mask, + static void apply(ArrayViewTuple&& arrayViews, const Mask& mask, const Function& function, const std::tuple& slicerArgs, const std::tuple& maskArgs) { // Iterate over this dimension. - if constexpr(Dim == ItrDim) { - + if constexpr (Dim == ItrDim) { // Get size of iteration dimenion from first view argument. const auto idxMax = std::get<0>(arrayViews).shape(ItrDim); forEach(idxMax, [&](idx_t idx) { - // Demote parallel execution policy to a non-parallel one in further recursion - ArrayForEachImpl, Dim + 1, ItrDims...>::apply( - std::forward(arrayViews), mask, function, - tuplePushBack(slicerArgs, idx), - tuplePushBack(maskArgs, idx)); + // Demote parallel execution policy to a non-parallel one in further + // recursion + ArrayForEachImpl< + execution::demote_policy_t, Dim + 1, + ItrDims...>::apply(std::forward(arrayViews), mask, + function, tuplePushBack(slicerArgs, idx), + tuplePushBack(maskArgs, idx)); }); } // Add a RangeAll to arguments. else { ArrayForEachImpl::apply( std::forward(arrayViews), mask, function, - tuplePushBack(slicerArgs, Range::all()), - maskArgs); + tuplePushBack(slicerArgs, Range::all()), maskArgs); } } }; template - struct is_applicable : std::false_type {}; +struct is_applicable : std::false_type {}; template -struct is_applicable> : std::is_invocable {}; +struct is_applicable> + : std::is_invocable {}; template -inline constexpr bool is_applicable_v = is_applicable::value; +inline constexpr bool is_applicable_v = is_applicable::value; template struct ArrayForEachImpl { - template - static void apply(ArrayViewTuple&& arrayViews, - const Mask& mask, + static void apply(ArrayViewTuple&& arrayViews, const Mask& mask, const Function& function, const std::tuple& slicerArgs, const std::tuple& maskArgs) { - constexpr auto maskPresent = !std::is_same_v; if constexpr (maskPresent) { - - constexpr auto invocableMask = std::is_invocable_r_v; - static_assert (invocableMask, - "Cannot invoke mask function with given arguments.\n" - "Make sure you arguments are N integers (or auto...) " - "where N == sizeof...(ItrDims). Function must return an int." - ); - - if (std::apply(mask, maskArgs)) { - return; - } - + constexpr auto invocableMask = + std::is_invocable_r_v; + static_assert( + invocableMask, + "Cannot invoke mask function with given arguments.\n" + "Make sure you arguments are N integers (or auto...) " + "where N == sizeof...(ItrDims). Function must return an int."); + + if (std::apply(mask, maskArgs)) { + return; + } } - auto slices = makeSlices(slicerArgs, std::forward(arrayViews)); + auto slices = + makeSlices(slicerArgs, std::forward(arrayViews)); - constexpr auto applicable = is_applicable_v; - static_assert(applicable, "Cannot invoke function with given arguments. " - "Make sure you the arguments are rvalue references (Slice&&) or const references (const Slice&) or regular value (Slice)" ); + constexpr auto applicable = is_applicable_v; + static_assert( + applicable, + "Cannot invoke function with given arguments. " + "Make sure you the arguments are rvalue references (Slice&&) or const " + "references (const Slice&) or regular value (Slice)"); std::apply(function, std::move(slices)); } - }; } // namespace detail @@ -286,41 +286,36 @@ struct ArrayForEach { /// and is executed with signature g(idx_i, idx_j,...), where the idxs /// are indices of ItrDims. /// When a config is supplied containing "execution_policy" = - /// "sequenced_policy" (default). All loops are then executed in sequential - /// (row-major) order. - /// With "execution_policy" = "parallel_unsequenced" the first loop is executed - /// using OpenMP. The remaining loops are executed in serial. - /// Note: The lowest ArrayView.rank() must be greater than or equal - /// to the highest dim in ItrDims. TODO: static checking for this. + /// "sequenced_policy" (default). All loops are then executed in + /// sequential (row-major) order. With "execution_policy" = + /// "parallel_unsequenced" the first loop is executed using OpenMP. + /// The remaining loops are executed in serial. Note: The lowest + /// ArrayView.rank() must be greater than or equal to the highest dim + /// in ItrDims. TODO: static checking for this. template static void apply(const eckit::Parametrisation& conf, - std::tuple&& arrayViews, - const Mask& mask, const Function& function) { - + std::tuple&& arrayViews, const Mask& mask, + const Function& function) { auto execute = [&](auto execution_policy) { apply(execution_policy, std::move(arrayViews), mask, function); }; using namespace execution; std::string execution_policy; - if (conf.get("execution_policy",execution_policy)) { + if (conf.get("execution_policy", execution_policy)) { if (execution_policy == policy_name(par_unseq)) { execute(par_unseq); - } - else if (execution_policy == policy_name(par)) { + } else if (execution_policy == policy_name(par)) { execute(par); - } - else if (execution_policy == policy_name(unseq)) { + } else if (execution_policy == policy_name(unseq)) { execute(unseq); - } - else if (execution_policy == policy_name(seq)) { + } else if (execution_policy == policy_name(seq)) { execute(seq); + } else { + throw_Exception("Unrecognized execution policy " + execution_policy, + Here()); } - else { - throw_Exception("Unrecognized execution policy "+execution_policy, Here()); - } - } - else { + } else { execute(seq); } } @@ -328,35 +323,47 @@ struct ArrayForEach { /// brief Apply "For-Each" method. /// /// details As above, but Execution policy is determined at compile-time. - template ()>> - static void apply(ExecutionPolicy, std::tuple&& arrayViews, const Mask& mask, const Function& function) { - detail::ArrayForEachImpl::apply( - std::move(arrayViews), mask, function, std::make_tuple(), std::make_tuple()); + template ()>> + static void apply(ExecutionPolicy, std::tuple&& arrayViews, + const Mask& mask, const Function& function) { + detail::ArrayForEachImpl::apply( + std::move(arrayViews), mask, function, std::make_tuple(), + std::make_tuple()); } /// brief Apply "For-Each" method /// /// details Apply ForEach with default execution policy. template - static void apply(std::tuple&& arrayViews, const Mask& mask, const Function& function) { - apply(std::move(arrayViews), mask, function); + static void apply(std::tuple&& arrayViews, const Mask& mask, + const Function& function) { + apply(std::move(arrayViews), mask, function); } /// brief Apply "For-Each" method /// - /// details Apply ForEach with run-time determined execution policy and no mask. + /// details Apply ForEach with run-time determined execution policy and no + /// mask. template - static void apply(const eckit::Parametrisation& conf, std::tuple&& arrayViews, const Function& function) { + static void apply(const eckit::Parametrisation& conf, + std::tuple&& arrayViews, + const Function& function) { apply(conf, std::move(arrayViews), detail::no_mask, function); } /// brief Apply "For-Each" method /// - /// details Apply ForEach with compile-time determined execution policy and no mask. + /// details Apply ForEach with compile-time determined execution policy and no + /// mask. template ()>> - static void apply(ExecutionPolicy executionPolicy, std::tuple&& arrayViews, const Function& function) { + typename = std::enable_if_t< + execution::is_execution_policy()>> + static void apply(ExecutionPolicy executionPolicy, + std::tuple&& arrayViews, + const Function& function) { apply(executionPolicy, std::move(arrayViews), detail::no_mask, function); } @@ -364,12 +371,22 @@ struct ArrayForEach { /// /// details Apply ForEach with default execution policy and no mask. template - static void apply(std::tuple&& arrayViews, const Function& function) { + static void apply(std::tuple&& arrayViews, + const Function& function) { apply(execution::seq, std::move(arrayViews), function); } - }; +/// brief Construct ArrayForEach and call apply +/// +/// details Construct an ArrayForEach using std::integer_sequence +/// . Remaining arguments are forwarded to apply +/// method. +template +void arrayForEachDim(std::integer_sequence, Args&&... args) { + ArrayForEach::apply(std::forward(args)...); +} + } // namespace helpers } // namespace array } // namespace atlas diff --git a/src/tests/array/test_array_foreach.cc b/src/tests/array/test_array_foreach.cc index f8c5b0b4a..42b9d2b09 100644 --- a/src/tests/array/test_array_foreach.cc +++ b/src/tests/array/test_array_foreach.cc @@ -1,5 +1,5 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * 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. @@ -7,13 +7,13 @@ #include #include +#include #include "atlas/array.h" #include "atlas/array/MakeView.h" #include "atlas/array/helpers/ArrayForEach.h" #include "atlas/array/helpers/ArraySlicer.h" #include "atlas/util/Config.h" - #include "tests/AtlasTestEnvironment.h" using namespace atlas::array; @@ -207,6 +207,61 @@ CASE("test_array_foreach_3_views") { EXPECT_EQ(count, 60); } +CASE("test_array_foreach_integer_sequence") { + + const auto arr1 = ArrayT(2, 3); + const auto view1 = make_view(arr1); + + const auto arr2 = ArrayT(2, 3, 4); + const auto view2 = make_view(arr2); + + const auto arr3 = ArrayT(2, 3, 4, 5); + const auto view3 = make_view(arr3); + + const auto zero = std::integer_sequence{}; + const auto one = std::integer_sequence{}; + const auto zeroOneTwoThree = std::make_integer_sequence{}; + + + // Test slice shapes. + + const auto loopFunctorDim0 = [](auto&& slice1, auto&& slice2, auto&& slice3) { + EXPECT_EQ(slice1.rank(), 1); + EXPECT_EQ(slice1.shape(0), 3); + + EXPECT_EQ(slice2.rank(), 2); + EXPECT_EQ(slice2.shape(0), 3); + EXPECT_EQ(slice2.shape(1), 4); + + EXPECT_EQ(slice3.rank(), 3); + EXPECT_EQ(slice3.shape(0), 3); + EXPECT_EQ(slice3.shape(1), 4); + EXPECT_EQ(slice3.shape(2), 5); + }; + arrayForEachDim(zero, std::tie(view1, view2, view3), loopFunctorDim0); + + const auto loopFunctorDim1 = [](auto&& slice1, auto&& slice2, auto&& slice3) { + EXPECT_EQ(slice1.rank(), 1); + EXPECT_EQ(slice1.shape(0), 2); + + EXPECT_EQ(slice2.rank(), 2); + EXPECT_EQ(slice2.shape(0), 2); + EXPECT_EQ(slice2.shape(1), 4); + + EXPECT_EQ(slice3.rank(), 3); + EXPECT_EQ(slice3.shape(0), 2); + EXPECT_EQ(slice3.shape(1), 4); + EXPECT_EQ(slice3.shape(2), 5); + }; + arrayForEachDim(one, std::tie(view1, view2, view3), loopFunctorDim1); + + // Test that slice resolves to double. + + const auto loopFunctorDimAll = [](auto&& slice3) { + static_assert(std::is_convertible_v); + }; + arrayForEachDim(zeroOneTwoThree, std::tie(view3), loopFunctorDimAll); +} CASE("test_array_foreach_forwarding") { From 8780370562c89ea81cf302d7aba72476844759c8 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 25 Jan 2024 16:42:27 +0000 Subject: [PATCH 08/23] Added test with empty integer sequence. --- src/tests/array/test_array_foreach.cc | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/tests/array/test_array_foreach.cc b/src/tests/array/test_array_foreach.cc index 42b9d2b09..303da5dd6 100644 --- a/src/tests/array/test_array_foreach.cc +++ b/src/tests/array/test_array_foreach.cc @@ -218,6 +218,7 @@ CASE("test_array_foreach_integer_sequence") { const auto arr3 = ArrayT(2, 3, 4, 5); const auto view3 = make_view(arr3); + const auto none = std::integer_sequence{}; const auto zero = std::integer_sequence{}; const auto one = std::integer_sequence{}; const auto zeroOneTwoThree = std::make_integer_sequence{}; @@ -225,6 +226,26 @@ CASE("test_array_foreach_integer_sequence") { // Test slice shapes. + const auto loopFunctorDimNone = [](auto&& slice1, auto&& slice2, + auto&& slice3) { + EXPECT_EQ(slice1.rank(), 2); + EXPECT_EQ(slice1.shape(0), 2); + EXPECT_EQ(slice1.shape(1), 3); + + EXPECT_EQ(slice2.rank(), 3); + EXPECT_EQ(slice2.shape(0), 2); + EXPECT_EQ(slice2.shape(1), 3); + EXPECT_EQ(slice2.shape(2), 4); + + EXPECT_EQ(slice3.rank(), 4); + EXPECT_EQ(slice3.shape(0), 2); + EXPECT_EQ(slice3.shape(1), 3); + EXPECT_EQ(slice3.shape(2), 4); + EXPECT_EQ(slice3.shape(3), 5); + }; + // No iterations. Apply functor directly to unsliced ArrayViews. + arrayForEachDim(none, std::tie(view1, view2, view3), loopFunctorDimNone); + const auto loopFunctorDim0 = [](auto&& slice1, auto&& slice2, auto&& slice3) { EXPECT_EQ(slice1.rank(), 1); EXPECT_EQ(slice1.shape(0), 3); From 61b2933258ce85a6ee41e9c74475f64fd6c61d30 Mon Sep 17 00:00:00 2001 From: Oliver Lomax Date: Tue, 27 Feb 2024 12:45:56 +0000 Subject: [PATCH 09/23] Refactoring of interpolation::method::SphericalVector and implementation of adjoint methods. (#168) * Fixed test output file name. * Added correct constness to geometry::UnitSphere::greatCircleCourse. * Removed unecessary RealMatrixMap definition. * Added arrayForEachDim wrapper function. * Partially implemented Complex Matrix multiply class. * Implimented ComplexMatrixMultiply class. * Further refactoring. * Added adjoint tests. * Adjusted test assert. * Attempting to solve MacOS nan issues. * Added norm checks on fields. * Assigning targetField = 0. in test. * Static cast in dot product. * Added explicit NaN checking. * Explicit types in NaN checking lambda parameters. * Checking for IEEE 754 arithmetic support. * Added NaN checking to weight generation. * Trying ArrayView assign instead of Array copy. * Disabled halo-exchange. * Addressed undefined behaviour with std::polar. * Fixed StructuredColumns inaccuracies. * Update SphericalVector.h Removed unnecessary headers. * Update SphericalVector.cc Removed unnecessary headers. * Update test_interpolation_spherical_vector.cc Removed initial error field assignment. * Added missing header include. * Updated copyright. * Added low-to-high resolution structured interpolation test. * Added high-to-low resolution structured interpolation test. * Restored test_array_foreach.cc. * Update ComplexMatrixMultiply.h Fixed strange formatting. * Addressed compilation failure when Eigen is disabled. * Added eckit Eigen macros to SparseMatrix.h. * Removed trial-and-error macros. * Implemented Eigen SparseMatrix iterators. * Separated out 2-vector and 3-vector MatMul methods. Replaced shared pointers with unique pointers. * Restored omp_parallel_for. * Updated self-destruct matrix class implemented when Eigen is disabled. * Forced consistency between integral types. * Replaced adjoint move operation with copy. * Added explicit on SparseMatrix private constructor. * Updated doxygen comment. SparseMatrix class no longe resembles eckit::SparseMatrix. * Removed unecessary copy * Addressed partial override warning. * Changed smart pointers to raw pointers. * Attempted to disable vectorisation if __NVCOMPILER is defined. * Fixed error on SparseMatrix dummy class. * Turned off OpenMP when __NVCOMPILER is defined. * Added stronger debug checks on weight matrices. Forced checks when __NVCOMPILER is defined. * I forgot to actually disable OpenMP where it mattered. I've disabled it properly now. * Removed unnecessary macros. * Changed ATLAS_ASSERT_MSG to ATLAS_ASSERT. * Fixed incorrect check for successful cast in SphericalVector constructor. * Removed SFINAE. * Replaced lambda expressions with private class methods. * Removed NVCOMPILER macros. * Removed some duplicated code in multiplyAdd method. --- src/atlas/CMakeLists.txt | 13 +- .../sphericalvector/ComplexMatrixMultiply.h | 217 ++++++++++++ .../method/sphericalvector/SparseMatrix.h | 107 ++++++ .../method/sphericalvector/SphericalVector.cc | 317 +++++++----------- .../method/sphericalvector/SphericalVector.h | 78 ++--- .../method/sphericalvector/Types.h | 32 ++ src/atlas/util/Geometry.h | 8 +- .../test_interpolation_spherical_vector.cc | 244 ++++++++++---- 8 files changed, 704 insertions(+), 312 deletions(-) create mode 100644 src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h create mode 100644 src/atlas/interpolation/method/sphericalvector/SparseMatrix.h create mode 100644 src/atlas/interpolation/method/sphericalvector/Types.h diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 3636d74d6..634a22dce 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -524,7 +524,7 @@ functionspace/detail/PointCloudInterface.cc functionspace/detail/CubedSphereStructure.h functionspace/detail/CubedSphereStructure.cc -# for cubedsphere matching mesh partitioner +# for cubedsphere matching mesh partitioner interpolation/method/cubedsphere/CellFinder.cc interpolation/method/cubedsphere/CellFinder.h interpolation/Vector2D.cc @@ -541,8 +541,8 @@ interpolation/element/Triag3D.cc interpolation/element/Triag3D.h interpolation/method/Intersect.cc interpolation/method/Intersect.h -interpolation/method/Ray.cc # For testing Quad -interpolation/method/Ray.h # For testing Quad +interpolation/method/Ray.cc # For testing Quad +interpolation/method/Ray.h # For testing Quad # for BuildConvexHull3D @@ -634,8 +634,11 @@ interpolation/method/knn/KNearestNeighboursBase.cc interpolation/method/knn/KNearestNeighboursBase.h interpolation/method/knn/NearestNeighbour.cc interpolation/method/knn/NearestNeighbour.h -interpolation/method/sphericalvector/SphericalVector.h +interpolation/method/sphericalvector/ComplexMatrixMultiply.h +interpolation/method/sphericalvector/SparseMatrix.h interpolation/method/sphericalvector/SphericalVector.cc +interpolation/method/sphericalvector/SphericalVector.h +interpolation/method/sphericalvector/Types.h interpolation/method/structured/Cubic2D.cc interpolation/method/structured/Cubic2D.h interpolation/method/structured/Cubic3D.cc @@ -868,7 +871,7 @@ if( NOT atlas_HAVE_ATLAS_FUNCTIONSPACE ) unset( atlas_parallel_srcs ) unset( atlas_output_srcs ) unset( atlas_redistribution_srcs ) - unset( atlas_linalg_srcs ) # only depends on array + unset( atlas_linalg_srcs ) # only depends on array endif() if( NOT atlas_HAVE_ATLAS_INTERPOLATION ) diff --git a/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h b/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h new file mode 100644 index 000000000..34cfdacb2 --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h @@ -0,0 +1,217 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * 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. + */ + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "atlas/array/ArrayView.h" +#include "atlas/array/LocalView.h" +#include "atlas/array/Range.h" +#include "atlas/array/helpers/ArrayForEach.h" +#include "atlas/interpolation/method/sphericalvector/Types.h" +#include "atlas/parallel/omp/omp.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace detail { + +struct TwoVectorTag {}; +struct ThreeVectorTag {}; +constexpr auto twoVector = TwoVectorTag{}; +constexpr auto threeVector = ThreeVectorTag{}; + +// Permits the equivalent of static_assert(false, msg). Will be addressed in +// C++26. +template +constexpr bool always_false_v = false; + +/// @brief Helper class to perform complex matrix multiplications +/// +/// @details Performs matrix multiplication between fields of 2-vectors or +/// 3-vectors. Fields must have Rank >= 2. Here, the assumption is +/// that Dim = 0 is the horizontal dimension, and Dim = (Rank - 1) is +/// the vector element dimension. +template +class ComplexMatrixMultiply { + public: + ComplexMatrixMultiply() = default; + + /// @brief Construct object from sparse matrices. + /// + /// @details complexWeights is a SparseMatrix of weights. realWeights is a + /// SparseMatrix containing the magnitudes of the elements of + /// complexWeights. + ComplexMatrixMultiply(const ComplexMatrix& complexWeights, + const RealMatrix& realWeights) + : complexWeightsPtr_{&complexWeights}, realWeightsPtr_{&realWeights} { + if constexpr (ATLAS_BUILD_TYPE_DEBUG) { + ATLAS_ASSERT(complexWeightsPtr_->rows() == realWeightsPtr_->rows()); + ATLAS_ASSERT(complexWeightsPtr_->cols() == realWeightsPtr_->cols()); + ATLAS_ASSERT(complexWeightsPtr_->nonZeros() == + realWeightsPtr_->nonZeros()); + for (auto rowIndex = Index{0}; rowIndex < complexWeightsPtr_->rows(); + ++rowIndex) { + for (auto [complexRowIter, realRowIter] = rowIters(rowIndex); + complexRowIter; ++complexRowIter, ++realRowIter) { + ATLAS_ASSERT(realRowIter); + ATLAS_ASSERT(realRowIter.row() == complexRowIter.row()); + ATLAS_ASSERT(realRowIter.col() == complexRowIter.col()); + // tinyNum ~= 2.3e-13 for double. + constexpr auto tinyNum = 1024 * std::numeric_limits::epsilon(); + const auto complexMagnitude = std::abs(complexRowIter.value()); + const auto realValue = realRowIter.value(); + const auto error = std::abs(complexMagnitude - realValue); + + const auto printError = [&]() { + auto strStream = std::ostringstream{}; + strStream << "Error complex components: "; + strStream << std::setprecision(19) << error; + return strStream.str(); + }; + ATLAS_ASSERT(error < tinyNum, printError()); + } + } + } + } + + /// @brief Apply complex matrix vector multiplication. + /// + /// @details Multiply weights by the elements in sourceView to give + /// elements in targetView. If VectorType == TwoVectorTag, + /// complexWeights are applied to the horizontal elements of + /// sourceView. If VectorType == ThreeVectorTag, then realWeights + /// are additionally applied to the vertical elements of sourceView. + template + void apply(const array::ArrayView& sourceView, + array::ArrayView& targetView, VectorType) const { + if constexpr (std::is_same_v) { + applyTwoVector(sourceView, targetView); + } else if constexpr (std::is_same_v) { + applyThreeVector(sourceView, targetView); + } else { + static_assert(always_false_v, "Unknown vector type"); + } + } + + private: + /// @brief Apply 2-vector MatMul. + template + void applyTwoVector(const array::ArrayView& sourceView, + array::ArrayView& targetView) const { + // We could probably optimise contiguous arrays using + // reinterpret_cast*>(view.data()). This is fine + // according to the C++ standard! + atlas_omp_parallel_for(auto rowIndex = Index{0}; + rowIndex < complexWeightsPtr_->rows(); ++rowIndex) { + auto targetSlice = sliceColumn(targetView, rowIndex); + if constexpr (InitialiseTarget) { + targetSlice.assign(0.); + } + + for (auto complexRowIter = complexWeightsPtr_->rowIter(rowIndex); + complexRowIter; ++complexRowIter) { + const auto colIndex = complexRowIter.col(); + const auto complexWeight = complexRowIter.value(); + const auto sourceSlice = sliceColumn(sourceView, colIndex); + multiplyAdd(sourceSlice, targetSlice, complexWeight); + } + } + } + + /// @brief Apply 3-vector MatMul. + template + void applyThreeVector(const array::ArrayView& sourceView, + array::ArrayView& targetView) const { + atlas_omp_parallel_for(auto rowIndex = Index{0}; + rowIndex < complexWeightsPtr_->rows(); ++rowIndex) { + auto targetSlice = sliceColumn(targetView, rowIndex); + if constexpr (InitialiseTarget) { + targetSlice.assign(0.); + } + + for (auto [complexRowIter, realRowIter] = rowIters(rowIndex); + complexRowIter; ++complexRowIter, ++realRowIter) { + const auto colIndex = complexRowIter.col(); + const auto complexWeight = complexRowIter.value(); + const auto realWeight = realRowIter.value(); + const auto sourceSlice = sliceColumn(sourceView, colIndex); + multiplyAdd(sourceSlice, targetSlice, complexWeight, realWeight); + } + } + } + + /// @brief Multiply source column by weight(s) and add to target column. + template + void multiplyAdd(const array::LocalView& sourceColumn, + array::LocalView& targetColumn, + Weights... weights) const { + for (auto levelIdx = 0; levelIdx < sourceColumn.shape(0); ++levelIdx) { + const auto sourceElem = sliceColumn(sourceColumn, levelIdx); + auto targetElem = sliceColumn(targetColumn, levelIdx); + multiplyAdd(sourceElem, targetElem, weights...); + } + } + + /// @brief Multiply source element by complex weight and add to target + /// element. + template + void multiplyAdd(const array::LocalView& sourceElem, + array::LocalView& targetElem, + Complex complexWeight) const { + const auto targetVector = + complexWeight * Complex(sourceElem(0), sourceElem(1)); + targetElem(0) += targetVector.real(); + targetElem(1) += targetVector.imag(); + } + + /// @brief Multiply source element by complex and real weights and add to + /// target element. + template + void multiplyAdd(const array::LocalView& sourceElem, + array::LocalView& targetElem, + Complex complexWeight, Real realWeight) const { + multiplyAdd(sourceElem, targetElem, complexWeight); + targetElem(2) += realWeight * sourceElem(2); + } + + /// @brief Return a pair of complex and real row iterators + std::pair rowIters( + Index rowIndex) const { + return std::make_pair(complexWeightsPtr_->rowIter(rowIndex), + realWeightsPtr_->rowIter(rowIndex)); + } + + /// @brief Makes the slice arrayView.slice(index, Range::all()...). + template + static auto sliceColumn(View& arrayView, Index index) { + constexpr auto Rank = std::decay_t::rank(); + using RangeAll = decltype(array::Range::all()); + + const auto slicerArgs = std::tuple_cat(std::make_tuple(index), + std::array{}); + const auto slicer = [&](auto&&... args) { + return arrayView.slice(args...); + }; + + return std::apply(slicer, slicerArgs); + } + + const ComplexMatrix* complexWeightsPtr_{}; + const RealMatrix* realWeightsPtr_{}; +}; + +} // namespace detail +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SparseMatrix.h b/src/atlas/interpolation/method/sphericalvector/SparseMatrix.h new file mode 100644 index 000000000..a399d80ca --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/SparseMatrix.h @@ -0,0 +1,107 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * 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. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/library/defines.h" +#if ATLAS_HAVE_EIGEN + +#include +#endif + +#include "atlas/runtime/Exception.h" +#include "eckit/log/CodeLocation.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace detail { + +#if ATLAS_HAVE_EIGEN +/// @brief Wrapper class for Eigen sparse matrix +/// +/// @details Eigen sparse matrix wrapper. Allows preprocessor disabling of class +/// if Eigen library is not present. +template +class SparseMatrix { + public: + using Index = int; + using EigenMatrix = Eigen::SparseMatrix; + using Triplet = Eigen::Triplet; + using Triplets = std::vector; + using RowIter = typename EigenMatrix::InnerIterator; + + SparseMatrix() = default; + + SparseMatrix(Index nRows, Index nCols, const Triplets& triplets) + : eigenMatrix_(nRows, nCols) { + eigenMatrix_.setFromTriplets(triplets.begin(), triplets.end()); + } + + Index nonZeros() const { return eigenMatrix_.nonZeros(); } + Index rows() const { return eigenMatrix_.rows(); } + Index cols() const { return eigenMatrix_.cols(); } + RowIter rowIter(Index rowIndex) const { + return RowIter(eigenMatrix_, rowIndex); + } + SparseMatrix adjoint() const { + auto adjointMatrix = SparseMatrix{}; + adjointMatrix.eigenMatrix_ = eigenMatrix_.adjoint(); + return adjointMatrix; + } + + private: + EigenMatrix eigenMatrix_{}; +}; +#else + +template +class SparseMatrix { + public: + using Index = int; + + class Triplet { + public: + template + constexpr Triplet(const Args&... args) {} + }; + using Triplets = std::vector; + + class RowIter { + public: + template + constexpr RowIter(const Args&... args) {} + constexpr Index row() const { return Index{}; } + constexpr Index col() const { return Index{}; } + constexpr Value value() const { return Value{}; } + constexpr operator bool() const { return false; } + constexpr RowIter& operator++() { return *this; } + }; + + constexpr SparseMatrix() = default; + template + SparseMatrix(const Args&... args) { + throw_Exception("Atlas has been compiled without Eigen", Here()); + } + constexpr Index nonZeros() const { return Index{}; } + constexpr Index rows() const { return Index{}; } + constexpr Index cols() const { return Index{}; } + constexpr RowIter rowIter(Index rowIndex) const { return RowIter{}; } + constexpr SparseMatrix adjoint() const { + return SparseMatrix{}; + } +}; +#endif + +} // namespace detail +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 483ce8634..4d46225ae 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -1,132 +1,48 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * 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. */ -#include "atlas/library/defines.h" #include "atlas/interpolation/method/sphericalvector/SphericalVector.h" -#include -#include -#include -#include - -#include "eckit/config/LocalConfiguration.h" - #include "atlas/array/ArrayView.h" -#include "atlas/array/helpers/ArrayForEach.h" -#include "atlas/array/Range.h" #include "atlas/field/Field.h" #include "atlas/field/FieldSet.h" #include "atlas/interpolation/Cache.h" #include "atlas/interpolation/Interpolation.h" #include "atlas/interpolation/method/MethodFactory.h" +#include "atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h" +#include "atlas/interpolation/method/sphericalvector/Types.h" +#include "atlas/option/Options.h" #include "atlas/parallel/omp/omp.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Trace.h" #include "atlas/util/Constants.h" #include "atlas/util/Geometry.h" +#include "eckit/config/LocalConfiguration.h" namespace atlas { namespace interpolation { namespace method { -namespace { -MethodBuilder __builder("spherical-vector"); -} - -#if ATLAS_HAVE_EIGEN - -// A bug exists in intel versions < intel/2022.2 with OpenMP -// Intel OneAPI version 2022.2 corresponds to Intel classic (icpc) version 2021.6 -#if defined(__INTEL_COMPILER) && defined(__INTEL_COMPILER_UPDATE) -#if (__INTEL_COMPILER <= 2021) && (__INTEL_COMPILER_UPDATE < 6) -#warning Disabling OpenMP to prevent internal compiler error for intel-classic version < 2021.6 (intel-oneapi/2022.2) -#undef atlas_omp_parallel_for -#define atlas_omp_parallel_for for -#endif -#endif - -using Complex = SphericalVector::Complex; - -template -using SparseMatrix = SphericalVector::SparseMatrix; -using RealMatrixMap = Eigen::Map>; -using ComplexTriplets = std::vector>; -using RealTriplets = std::vector>; - -using EckitMatrix = eckit::linalg::SparseMatrix; - namespace { - -RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { - return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(), - baseMatrix.nonZeros(), baseMatrix.outer(), - baseMatrix.inner(), baseMatrix.data()); -} - -template -auto getInnerIt(const Matrix& matrix, typename Matrix::Index k) { - return typename Matrix::InnerIterator(matrix, k); -} - -template -void sparseMatrixForEach(const Functor& functor, const Matrix& matrix) { - using Index = typename Matrix::Index; - atlas_omp_parallel_for (auto k = Index{}; k < matrix.outerSize(); ++k) { - for (auto it = getInnerIt(matrix, k); it; ++it) { - functor(it.row(), it.col(), it.value()); - } - } +MethodBuilder __builder("spherical-vector"); } -template -void sparseMatrixForEach(const Functor& functor, const Matrix1& matrix1, - const Matrix2& matrix2) { - using Index = typename Matrix1::Index; - atlas_omp_parallel_for (auto k = Index{}; k < matrix1.outerSize(); ++k) { - for (auto [it1, it2] = - std::make_pair(getInnerIt(matrix1, k), getInnerIt(matrix2, k)); - it1; ++it1, ++it2) { - functor(it1.row(), it1.col(), it1.value(), it2.value()); - } - } -} +using namespace detail; -template -void matrixMultiply(const SourceView& sourceView, TargetView& targetView, - const Functor& multiplyFunctor, - const Matrices&... matrices) { - - const auto multiplyColumn = [&](auto i, auto j, const auto&... weights) { - constexpr auto Rank = std::decay_t::rank(); - if constexpr (Rank == 2) { - const auto sourceSlice = sourceView.slice(j, array::Range::all()); - auto targetSlice = targetView.slice(i, array::Range::all()); - multiplyFunctor(sourceSlice, targetSlice, weights...); - } else if constexpr(Rank == 3) { - const auto sourceSlice = - sourceView.slice(j, array::Range::all(), array::Range::all()); - auto targetSlice = - targetView.slice(i, array::Range::all(), array::Range::all()); - array::helpers::ArrayForEach<0>::apply( - std::tie(sourceSlice, targetSlice), - [&](auto&& sourceVars, auto&& targetVars) { - multiplyFunctor(sourceVars, targetVars, weights...); - }); - } else { - ATLAS_NOTIMPLEMENTED; - } - }; +using WeightsMatMul = detail::ComplexMatrixMultiply; +using WeightsMatMulAdjoint = detail::ComplexMatrixMultiply; - sparseMatrixForEach(multiplyColumn, matrices...); +SphericalVector::SphericalVector(const Config& config) : Method(config) { + const auto* conf = dynamic_cast(&config); + ATLAS_ASSERT(conf, "config must be derived from eckit::LocalConfiguration"); + interpolationScheme_ = conf->getSubConfiguration("scheme"); + adjoint_ = conf->getBool("adjoint", false); } -} // namespace - void SphericalVector::do_setup(const Grid& source, const Grid& target, const Cache&) { ATLAS_NOTIMPLEMENTED; @@ -145,52 +61,66 @@ void SphericalVector::do_setup(const FunctionSpace& source, setMatrix(Interpolation(interpolationScheme_, source_, target_)); // Get matrix data. - const auto nRows = matrix().rows(); - const auto nCols = matrix().cols(); - const auto nNonZeros = matrix().nonZeros(); - const auto baseWeights = makeMatrixMap(matrix()); + const auto nRows = static_cast(matrix().rows()); + const auto nCols = static_cast(matrix().cols()); + const auto nNonZeros = static_cast(matrix().nonZeros()); + const auto* outerIndices = matrix().outer(); + const auto* innerIndices = matrix().inner(); + const auto* baseWeights = matrix().data(); // Note: need to store copy of weights as Eigen3 sorts compressed rows by j // whereas eckit does not. - complexWeights_ = std::make_shared(nRows, nCols); - realWeights_ = std::make_shared(nRows, nCols); auto complexTriplets = ComplexTriplets(nNonZeros); auto realTriplets = RealTriplets(nNonZeros); - const auto sourceLonLats = array::make_view(source_.lonlat()); - const auto targetLonLats = array::make_view(target_.lonlat()); - - geometry::UnitSphere unitSphere; - - const auto setWeights = [&](auto i, auto j, const auto& baseWeight) { - const auto sourceLonLat = - PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); - const auto targetLonLat = - PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); - - const auto alpha = unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); - - const auto deltaAlpha = - (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - - const auto idx = std::distance(baseWeights.valuePtr(), &baseWeight); - - complexTriplets[idx] = {int(i), int(j), std::polar(baseWeight, deltaAlpha)}; - realTriplets[idx] = {int(i), int(j), baseWeight}; - }; - - sparseMatrixForEach(setWeights, baseWeights); - complexWeights_->setFromTriplets(complexTriplets.begin(), - complexTriplets.end()); - realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end()); + // Make sure halo lonlats are same as owned points. + auto sourceLonLats = source_.createField(option::name("lonlat") | + option::variables(2)); + auto targetLonLats = target_.createField(option::name("lonlat") | + option::variables(2)); + sourceLonLats.array().copy(source_.lonlat()); + targetLonLats.array().copy(target_.lonlat()); + sourceLonLats.haloExchange(); + targetLonLats.haloExchange(); + + const auto sourceLonLatsView = array::make_view(sourceLonLats); + const auto targetLonLatsView = array::make_view(targetLonLats); + + const auto unitSphere = geometry::UnitSphere{}; + + atlas_omp_parallel_for(auto rowIndex = Index{0}; rowIndex < nRows; + ++rowIndex) { + for (auto dataIndex = outerIndices[rowIndex]; + dataIndex < outerIndices[rowIndex + 1]; ++dataIndex) { + const auto colIndex = static_cast(innerIndices[dataIndex]); + const auto baseWeight = baseWeights[dataIndex]; + + const auto sourceLonLat = PointLonLat(sourceLonLatsView(colIndex, 0), + sourceLonLatsView(colIndex, 1)); + const auto targetLonLat = PointLonLat(targetLonLatsView(rowIndex, 0), + targetLonLatsView(rowIndex, 1)); + + const auto alpha = + unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); + + const auto deltaAlpha = + (alpha.first - alpha.second) * util::Constants::degreesToRadians(); + + complexTriplets[dataIndex] = + ComplexTriplet{rowIndex, colIndex, + Complex{baseWeight * std::cos(deltaAlpha), + baseWeight * std::sin(deltaAlpha)}}; + realTriplets[dataIndex] = RealTriplet{rowIndex, colIndex, baseWeight}; + } + } - ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); - ATLAS_ASSERT(realWeights_->nonZeros() == matrix().nonZeros()); -} + complexWeights_ = ComplexMatrix(nRows, nCols, complexTriplets); + realWeights_ = RealMatrix(nRows, nCols, realTriplets); -SphericalVector::SphericalVector(const Config& config) : Method(config) { - const auto& conf = dynamic_cast(config); - interpolationScheme_ = conf.getSubConfiguration("scheme"); + if (adjoint_) { + complexWeightsAdjoint_ = complexWeights_.adjoint(); + realWeightsAdjoint_ = realWeights_.adjoint(); + } } void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } @@ -199,9 +129,7 @@ void SphericalVector::do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, Metadata& metadata) const { ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); - - const auto nFields = sourceFieldSet.size(); - ATLAS_ASSERT(nFields == targetFieldSet.size()); + ATLAS_ASSERT(sourceFieldSet.size() == targetFieldSet.size()); for (auto i = 0; i < sourceFieldSet.size(); ++i) { do_execute(sourceFieldSet[i], targetFieldSet[i], metadata); @@ -211,107 +139,118 @@ void SphericalVector::do_execute(const FieldSet& sourceFieldSet, void SphericalVector::do_execute(const Field& sourceField, Field& targetField, Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); - const auto fieldType = sourceField.metadata().getString("type", ""); if (fieldType != "vector") { - auto metadata = Metadata(); Method::do_execute(sourceField, targetField, metadata); return; } - if (target_.size() == 0) { - return; - } - - ATLAS_ASSERT_MSG(sourceField.variables() == 2 || sourceField.variables() == 3, - "Vector field can only have 2 or 3 components."); - Method::check_compatibility(sourceField, targetField, matrix()); haloExchange(sourceField); - - if (sourceField.datatype().kind() == array::DataType::KIND_REAL64) { - interpolate_vector_field(sourceField, targetField); - } else if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { - interpolate_vector_field(sourceField, targetField); - } else { - ATLAS_NOTIMPLEMENTED; - } - + interpolate_vector_field(sourceField, targetField, + WeightsMatMul(complexWeights_, realWeights_)); targetField.set_dirty(); } void SphericalVector::do_execute_adjoint(FieldSet& sourceFieldSet, const FieldSet& targetFieldSet, Metadata& metadata) const { - ATLAS_NOTIMPLEMENTED; + ATLAS_TRACE( + "atlas::interpolation::method::SphericalVector::do_execute_adjoint()"); + ATLAS_ASSERT(sourceFieldSet.size() == targetFieldSet.size()); + + for (auto i = 0; i < sourceFieldSet.size(); ++i) { + do_execute_adjoint(sourceFieldSet[i], targetFieldSet[i], metadata); + } } void SphericalVector::do_execute_adjoint(Field& sourceField, const Field& targetField, Metadata& metadata) const { - ATLAS_NOTIMPLEMENTED; + ATLAS_TRACE( + "atlas::interpolation::method::SphericalVector::do_execute_adjoint()"); + + const auto fieldType = sourceField.metadata().getString("type", ""); + if (fieldType != "vector") { + auto metadata = Metadata(); + Method::do_execute_adjoint(sourceField, targetField, metadata); + + return; + } + + Method::check_compatibility(sourceField, targetField, matrix()); + + ATLAS_ASSERT(adjoint_, "\"adjoint\" needs to be set to \"true\" in Config."); + interpolate_vector_field( + targetField, sourceField, + WeightsMatMulAdjoint(complexWeightsAdjoint_, realWeightsAdjoint_)); + adjointHaloExchange(sourceField); } -template +template +void SphericalVector::interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul) { + if (targetField.size() == 0) { + return; + } + + ATLAS_ASSERT_MSG(sourceField.variables() == 2 || sourceField.variables() == 3, + "Vector field can only have 2 or 3 components."); + + if (sourceField.datatype().kind() == array::DataType::KIND_REAL64) { + interpolate_vector_field(sourceField, targetField, matMul); + return; + } + + if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { + interpolate_vector_field(sourceField, targetField, matMul); + return; + } + + ATLAS_NOTIMPLEMENTED; +}; + +template void SphericalVector::interpolate_vector_field(const Field& sourceField, - Field& targetField) const { + Field& targetField, + const MatMul& matMul) { if (sourceField.rank() == 2) { - interpolate_vector_field(sourceField, targetField); + interpolate_vector_field(sourceField, targetField, matMul); return; } if (sourceField.rank() == 3) { - interpolate_vector_field(sourceField, targetField); + interpolate_vector_field(sourceField, targetField, matMul); return; } ATLAS_NOTIMPLEMENTED; } -template +template void SphericalVector::interpolate_vector_field(const Field& sourceField, - Field& targetField) const { - + Field& targetField, + const MatMul& matMul) { const auto sourceView = array::make_view(sourceField); auto targetView = array::make_view(targetField); - targetView.assign(0.); - - const auto horizontalComponent = [](const auto& sourceVars, auto& targetVars, - const auto& complexWeight) { - const auto sourceVector = Complex(sourceVars(0), sourceVars(1)); - const auto targetVector = complexWeight * sourceVector; - targetVars(0) += targetVector.real(); - targetVars(1) += targetVector.imag(); - }; if (sourceField.variables() == 2) { - matrixMultiply(sourceView, targetView, horizontalComponent, - *complexWeights_); + matMul.apply(sourceView, targetView, twoVector); return; } if (sourceField.variables() == 3) { - - const auto horizontalAndVerticalComponent = [&]( - const auto& sourceVars, auto& targetVars, const auto& complexWeight, - const auto& realWeight) { - horizontalComponent(sourceVars, targetVars, complexWeight); - targetVars(2) += realWeight * sourceVars(2); - }; - - matrixMultiply(sourceView, targetView, horizontalAndVerticalComponent, - *complexWeights_, *realWeights_); + matMul.apply(sourceView, targetView, threeVector); return; } ATLAS_NOTIMPLEMENTED; } -#endif - } // namespace method } // namespace interpolation } // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index 91ce2a1a2..c3514b0d1 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -1,41 +1,22 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * 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. */ -#include "atlas/library/defines.h" - #pragma once -#include -#include - -#if ATLAS_HAVE_EIGEN -#include -#endif - #include "atlas/functionspace/FunctionSpace.h" #include "atlas/interpolation/method/Method.h" -#include "atlas/linalg/sparse.h" - +#include "atlas/interpolation/method/sphericalvector/Types.h" namespace atlas { namespace interpolation { namespace method { -#if ATLAS_HAVE_EIGEN class SphericalVector : public Method { public: - using Complex = std::complex; - - template - using SparseMatrix = Eigen::SparseMatrix; - using ComplexMatrix = SparseMatrix; - using RealMatrix = SparseMatrix; - - /// @brief Interpolation post-processor for vector field interpolation /// /// @details Takes a base interpolation config keyed to "scheme" and creates @@ -70,50 +51,37 @@ class SphericalVector : public Method { Metadata& metadata) const override; private: - template - void interpolate_vector_field(const Field& source, Field& target) const; + template + static void interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul); - template - void interpolate_vector_field(const Field& source, Field& target) const; + template + static void interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul); + + template + static void interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul); using Method::do_setup; void do_setup(const FunctionSpace& source, const FunctionSpace& target) override; void do_setup(const Grid& source, const Grid& target, const Cache&) override; - eckit::LocalConfiguration interpolationScheme_; - - FunctionSpace source_; - FunctionSpace target_; + eckit::LocalConfiguration interpolationScheme_{}; + bool adjoint_{}; - std::shared_ptr complexWeights_; - std::shared_ptr realWeights_; + FunctionSpace source_{}; + FunctionSpace target_{}; + detail::ComplexMatrix complexWeights_{}; + detail::RealMatrix realWeights_{}; + detail::ComplexMatrix complexWeightsAdjoint_{}; + detail::RealMatrix realWeightsAdjoint_{}; }; -#else - class SphericalVector : public Method { - public: - SphericalVector(const Config& config) : Method(config) { - ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen"); - } - - ~SphericalVector() override {} - - void print(std::ostream&) const override {} - const FunctionSpace& source() const override {ATLAS_NOTIMPLEMENTED;} - const FunctionSpace& target() const override {ATLAS_NOTIMPLEMENTED;} - - void do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, - Metadata& metadata) const override {} - void do_execute(const Field& sourceField, Field& targetField, - Metadata& metadata) const override {} - private: - void do_setup(const FunctionSpace& source, - const FunctionSpace& target) override {} - void do_setup(const Grid& source, const Grid& target, const Cache&) override {} - }; -#endif - } // namespace method } // namespace interpolation diff --git a/src/atlas/interpolation/method/sphericalvector/Types.h b/src/atlas/interpolation/method/sphericalvector/Types.h new file mode 100644 index 000000000..f83e9640b --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/Types.h @@ -0,0 +1,32 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * 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. + */ + +#pragma once + +#include + +#include "atlas/interpolation/method/sphericalvector/SparseMatrix.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace detail { + +using Real = double; +using Complex = std::complex; +using ComplexMatrix = SparseMatrix; +using RealMatrix = SparseMatrix; +using Index = ComplexMatrix::Index; +using ComplexTriplet = ComplexMatrix::Triplet; +using ComplexTriplets = ComplexMatrix::Triplets; +using RealTriplet = RealMatrix::Triplet; +using RealTriplets = RealMatrix::Triplets; + +} // detail +} // method +} // interpolation +} // atlas diff --git a/src/atlas/util/Geometry.h b/src/atlas/util/Geometry.h index dd8173c3a..90b868ce8 100644 --- a/src/atlas/util/Geometry.h +++ b/src/atlas/util/Geometry.h @@ -51,7 +51,7 @@ class GeometryBase : public util::Object { virtual double distance(const Point3& p1, const Point3& p2) const = 0; virtual double radius() const = 0; virtual double area() const = 0; - virtual std::pair greatCircleCourse(const Point2& p1, const Point2& p2) = 0; + virtual std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const = 0; Point3 xyz(const Point2& lonlat) const { Point3 xyz; @@ -80,7 +80,7 @@ class GeometrySphereT : public GeometryBase { double distance(const Point3& p1, const Point3& p2) const override { return SphereT::distance(p1, p2); } double radius() const override { return SphereT::radius(); } double area() const override { return SphereT::area(); } - std::pair greatCircleCourse(const Point2& p1, const Point2& p2) override { + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const override { return atlas::geometry::detail::greatCircleCourse(p1, p2); } }; @@ -96,7 +96,7 @@ class GeometrySphere : public GeometryBase { double distance(const Point3& p1, const Point3& p2) const override { return Sphere::distance(radius_, p1, p2); } double radius() const override { return radius_; } double area() const override { return Sphere::area(radius_); } - std::pair greatCircleCourse(const Point2& p1, const Point2& p2) override { + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const override { return atlas::geometry::detail::greatCircleCourse(p1, p2); } @@ -129,7 +129,7 @@ class Geometry : public util::ObjectHandle { double distance(const Point3& p1, const Point3& p2) const { return get()->distance(p1, p2); } double radius() const { return get()->radius(); } double area() const { return get()->area(); } - std::pair greatCircleCourse(const Point2& p1, const Point2& p2) { return get()->greatCircleCourse(p1, p2); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const { return get()->greatCircleCourse(p1, p2); } protected: template diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index 1e68b8e91..19c778fdd 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -1,10 +1,14 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * 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. */ +#include +#include +#include + #include "atlas/array.h" #include "atlas/array/helpers/ArrayForEach.h" #include "atlas/field.h" @@ -58,13 +62,17 @@ double vortexVertical(double lon, double lat) { void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { - const auto& functionSpace = fieldSet[0].functionspace(); - const auto& mesh = functionspace::NodeColumns(functionSpace).mesh(); + const auto functionSpace = fieldSet[0].functionspace(); + const auto structuredColums = functionspace::StructuredColumns(functionSpace); + const auto nodeColumns = functionspace::NodeColumns(functionSpace); + const auto mesh = + structuredColums ? Mesh(structuredColums.grid()) : nodeColumns.mesh(); - const auto gmshConfig = util::Config("coordinates", "xyz") | - util::Config("ghost", true) | - util::Config("info", true); + const auto gmshConfig = Config("coordinates", "xyz") | Config("ghost", true) | + Config("info", true); const auto gmsh = output::Gmsh(fileName, gmshConfig); + + gmsh.write(mesh); gmsh.write(fieldSet, functionSpace); } @@ -80,10 +88,17 @@ const auto generateNodeColums(const std::string& gridName, // Helper struct to key different Functionspaces to strings struct FunctionSpaceFixtures { static const FunctionSpace& get(const std::string& fixture) { - static std::map functionSpaces = { - {"cubedsphere_mesh", - generateNodeColums("CS-LFR-48", "cubedsphere_dual")}, - {"gaussian_mesh", generateNodeColums("O48", "structured")}}; + static const auto functionSpaces = + std::map{ + {"cubedsphere_mesh", + generateNodeColums("CS-LFR-48", "cubedsphere_dual")}, + {"gaussian_mesh", generateNodeColums("O48", "structured")}, + {"structured_columns", + functionspace::StructuredColumns(Grid("O48"), option::halo(1))}, + {"structured_columns_lowres", + functionspace::StructuredColumns(Grid("O24"), option::halo(1))}, + {"structured_columns_hires", + functionspace::StructuredColumns(Grid("O96"), option::halo(1))}}; return functionSpaces.at(fixture); } }; @@ -91,7 +106,7 @@ struct FunctionSpaceFixtures { // Helper struct to key different grid configs to strings struct FieldSpecFixtures { static const Config& get(const std::string& fixture) { - static std::map fieldSpecs = { + static const auto fieldSpecs = std::map{ {"2vector", option::name("test field") | option::variables(2) | option::type("vector")}, {"3vector", option::name("test field") | option::variables(3) | @@ -100,6 +115,54 @@ struct FieldSpecFixtures { } }; +// Helper stcut to key different interpolation schemes to strings +struct InterpSchemeFixtures { + static const Config& get(const std::string& fixture) { + + static const auto cubedsphereBilinear = + option::type("cubedsphere-bilinear"); + static const auto finiteElement = option::type("finite-element"); + static const auto structuredLinear = + option::type("structured-linear2D") | option::halo(1); + + static const auto sphericalVector = + option::type("spherical-vector") | Config("adjoint", true); + + static const auto interpSchemes = std::map{ + {"cubedsphere_bilinear", cubedsphereBilinear}, + {"finite_element", finiteElement}, + {"structured_linear", structuredLinear}, + {"cubedsphere_bilinear_spherical", + sphericalVector | Config("scheme", cubedsphereBilinear)}, + {"finite_element_spherical", + sphericalVector | Config("scheme", finiteElement)}, + {"structured_linear_spherical", + sphericalVector | Config("scheme", structuredLinear)}}; + return interpSchemes.at(fixture); + } +}; + +template +double dotProduct(const array::ArrayView& a, + const array::ArrayView& b) { + auto dotProd = 0.; + arrayForEachDim(std::make_integer_sequence{}, std::tie(a, b), + [&](const double& aElem, + const double& bElem) { dotProd += aElem * bElem; }); + return dotProd; +} + +template +int countNans(const array::ArrayView& view) { + auto nNans = 0; + arrayForEachDim(std::make_integer_sequence{}, std::tie(view), + [&](const double& viewElem) { + if (!std::isfinite(viewElem)) { + ++nNans; + } + }); + return nNans; +} template void testInterpolation(const Config& config) { @@ -121,13 +184,15 @@ void testInterpolation(const Config& config) { FieldSpecFixtures::get(config.getString("field_spec_fixture")); if constexpr (Rank == 3) fieldSpec.set("levels", 2); - auto sourceField = array::make_view( - sourceFieldSet.add(sourceFunctionSpace.createField(fieldSpec))); + auto sourceField = + sourceFieldSet.add(sourceFunctionSpace.createField(fieldSpec)); + auto targetField = + targetFieldSet.add(targetFunctionSpace.createField(fieldSpec)); - auto targetField = array::make_view( - targetFieldSet.add(targetFunctionSpace.createField(fieldSpec))); + auto sourceView = array::make_view(sourceField); + auto targetView = array::make_view(targetField); - ArrayForEach<0>::apply(std::tie(sourceLonLat, sourceField), + ArrayForEach<0>::apply(std::tie(sourceLonLat, sourceView), [](auto&& lonLat, auto&& sourceColumn) { const auto setElems = [&](auto&& sourceElem) { @@ -144,20 +209,22 @@ void testInterpolation(const Config& config) { }); sourceFieldSet.set_dirty(false); - const auto interp = Interpolation(config.getSubConfiguration("scheme"), - sourceFunctionSpace, targetFunctionSpace); + const auto interp = Interpolation( + InterpSchemeFixtures::get(config.getString("interp_fixture")), + sourceFunctionSpace, targetFunctionSpace); interp.execute(sourceFieldSet, targetFieldSet); targetFieldSet.haloExchange(); auto errorFieldSpec = fieldSpec; - errorFieldSpec.remove("variables"); + errorFieldSpec.remove("variables").set("name", "error field"); - auto errorField = array::make_view(targetFieldSet.add( + auto errorView = array::make_view(targetFieldSet.add( targetFunctionSpace.createField(errorFieldSpec))); + errorView.assign(0.); auto maxError = 0.; - ArrayForEach<0>::apply(std::tie(targetLonLat, targetField, errorField), + ArrayForEach<0>::apply(std::tie(targetLonLat, targetView, errorView), [&](auto&& lonLat, auto&& targetColumn, auto&& errorColumn) { @@ -179,62 +246,121 @@ void testInterpolation(const Config& config) { maxError = std::max(maxError, static_cast(errorElem)); }; - if constexpr(Rank == 2) { calcError(targetColumn, errorColumn); } - else if constexpr (Rank == 3) { + if + constexpr(Rank == 2) { calcError(targetColumn, errorColumn); } + else if + constexpr(Rank == 3) { ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), calcError); - } + } }); EXPECT_APPROX_EQ(maxError, 0., config.getDouble("tol")); gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); + + + // Adjoint test + auto targetAdjoint = targetFunctionSpace.createField(fieldSpec); + auto targetAdjointView = array::make_view(targetAdjoint); + targetAdjoint.array().copy(targetField); + targetAdjoint.adjointHaloExchange(); + + auto sourceAdjoint = sourceFunctionSpace.createField(fieldSpec); + auto sourceAdjointView = array::make_view(sourceAdjoint); + sourceAdjointView.assign(0.); + sourceAdjoint.set_dirty(false); + + interp.execute_adjoint(sourceAdjoint, targetAdjoint); + + // Check fields for nans or +/-inf + EXPECT_EQ(countNans(targetView), 0); + EXPECT_EQ(countNans(sourceView), 0); + EXPECT_EQ(countNans(targetAdjointView), 0); + EXPECT_EQ(countNans(sourceAdjointView), 0); + + constexpr auto tinyNum = 1e-13; + const auto targetDotTarget = dotProduct(targetView, targetView); + const auto sourceDotSourceAdjoint = dotProduct(sourceView, sourceAdjointView); + const auto dotProdRatio = targetDotTarget / sourceDotSourceAdjoint; + EXPECT_APPROX_EQ(dotProdRatio, 1., tinyNum); } CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { - - const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear"); - const auto interpScheme = - util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); - const auto cubedSphereConf = Config("source_fixture", "cubedsphere_mesh") - .set("target_fixture", "gaussian_mesh") - .set("field_spec_fixture", "2vector") - .set("file_id", "spherical_vector_cs2") - .set("scheme", interpScheme) - .set("tol", 0.00018); - - testInterpolation((cubedSphereConf)); + const auto config = + Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "cubedsphere_bilinear_spherical") + .set("file_id", "spherical_vector_cs2") + .set("tol", 0.00018); + + testInterpolation((config)); } CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { - - const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear"); - const auto interpScheme = - util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); - const auto cubedSphereConf = Config("source_fixture", "cubedsphere_mesh") - .set("target_fixture", "gaussian_mesh") - .set("field_spec_fixture", "3vector") - .set("file_id", "spherical_vector_cs3") - .set("scheme", interpScheme) - .set("tol", 0.00096); - - testInterpolation((cubedSphereConf)); + const auto config = + Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "3vector") + .set("interp_fixture", "cubedsphere_bilinear_spherical") + .set("file_id", "spherical_vector_cs3") + .set("tol", 0.00096); + + testInterpolation((config)); } CASE("finite element vector interpolation (2d-field, 2-vector)") { + const auto config = + Config("source_fixture", "gaussian_mesh") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "finite_element_spherical") + .set("file_id", "spherical_vector_fe") + .set("tol", 0.00015); + + testInterpolation((config)); +} + +CASE("structured columns vector interpolation (2d-field, 2-vector)") { + + const auto config = + Config("source_fixture", "structured_columns") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc") + .set("tol", 0.00017); - const auto baseInterpScheme = util::Config("type", "finite-element"); - const auto interpScheme = - util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); - const auto cubedSphereConf = Config("source_fixture", "gaussian_mesh") - .set("target_fixture", "cubedsphere_mesh") - .set("field_spec_fixture", "2vector") - .set("file_id", "spherical_vector_cs") - .set("scheme", interpScheme) - .set("tol", 0.00015); - - testInterpolation((cubedSphereConf)); + testInterpolation((config)); } + +CASE("structured columns vector interpolation (2d-field, 2-vector, low-res)") { + + const auto config = + Config("source_fixture", "structured_columns_lowres") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_lr") + .set("tol", 0.00056); + + testInterpolation((config)); +} + +CASE("structured columns vector interpolation (2d-field, 2-vector, hi-res)") { + + const auto config = + Config("source_fixture", "structured_columns_hires") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_hr") + .set("tol", 0.000044); + + testInterpolation((config)); +} + } } From 4e2a05546bd32ed6762da9fcf69ca4a1978dc72e Mon Sep 17 00:00:00 2001 From: Benjamin Menetrier Date: Tue, 27 Feb 2024 10:18:25 +0100 Subject: [PATCH 10/23] Bugfix for regional grids with ny > nx --- src/atlas/functionspace/detail/StructuredColumns_setup.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/functionspace/detail/StructuredColumns_setup.cc b/src/atlas/functionspace/detail/StructuredColumns_setup.cc index 01a9a1be4..4bdc0a4e9 100644 --- a/src/atlas/functionspace/detail/StructuredColumns_setup.cc +++ b/src/atlas/functionspace/detail/StructuredColumns_setup.cc @@ -371,7 +371,7 @@ void StructuredColumns::setup(const grid::Distribution& distribution, const ecki idx_t jj_max = j + halo; if (regional) { jj_min = std::max(jj_min, idx_t{0}); - jj_max = std::min(jj_max, grid_->nx(j) - idx_t{1}); + jj_max = std::min(jj_max, grid_->ny() - idx_t{1}); } for (idx_t jj = jj_min; jj <= jj_max; ++jj) { idx_t jjj = compute_j(jj); From c3333e96dba6f432ac3c83b8f69c75e578e57df6 Mon Sep 17 00:00:00 2001 From: Oliver Lomax Date: Wed, 6 Mar 2024 08:57:52 +0000 Subject: [PATCH 11/23] Update `SphericalVector` to work with StructuredColumns as source functionspace. (#175) * Realised StructuredColumns weight corrections cancel due to Euler identity. * Added great-circle course tests for points on the pole. * Added great-circle course special case where lonlat1 == lonlat2. --- .../method/sphericalvector/SphericalVector.cc | 9 +- src/atlas/util/Geometry.cc | 116 ++++++------ .../test_interpolation_spherical_vector.cc | 174 ++++++++---------- src/tests/util/test_unitsphere.cc | 28 ++- 4 files changed, 162 insertions(+), 165 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 4d46225ae..a8cc3a96f 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -78,13 +78,8 @@ void SphericalVector::do_setup(const FunctionSpace& source, option::variables(2)); auto targetLonLats = target_.createField(option::name("lonlat") | option::variables(2)); - sourceLonLats.array().copy(source_.lonlat()); - targetLonLats.array().copy(target_.lonlat()); - sourceLonLats.haloExchange(); - targetLonLats.haloExchange(); - - const auto sourceLonLatsView = array::make_view(sourceLonLats); - const auto targetLonLatsView = array::make_view(targetLonLats); + const auto sourceLonLatsView = array::make_view(source_.lonlat()); + const auto targetLonLatsView = array::make_view(target_.lonlat()); const auto unitSphere = geometry::UnitSphere{}; diff --git a/src/atlas/util/Geometry.cc b/src/atlas/util/Geometry.cc index 074d72d34..a274d2ead 100644 --- a/src/atlas/util/Geometry.cc +++ b/src/atlas/util/Geometry.cc @@ -3,7 +3,7 @@ * * 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 immGeometryies + * 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. */ @@ -12,9 +12,6 @@ #include -#include "eckit/geometry/Point2.h" -#include "eckit/geometry/Point3.h" - #include "atlas/library/config.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Constants.h" @@ -25,16 +22,15 @@ namespace geometry { namespace detail { void GeometrySphere::lonlat2xyz(const Point2& lonlat, Point3& xyz) const { #if ATLAS_ECKIT_VERSION_AT_LEAST(1, 24, 0) - Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true); + Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true); #else - Sphere::convertSphericalToCartesian(radius_, lonlat, xyz); + Sphere::convertSphericalToCartesian(radius_, lonlat, xyz); #endif } void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const { - Sphere::convertCartesianToSpherical(radius_, xyz, lonlat); + Sphere::convertCartesianToSpherical(radius_, xyz, lonlat); } - /// @brief Calculate great-cricle course between points /// /// @details Calculates the direction (clockwise from north) of a great-circle @@ -45,7 +41,10 @@ void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const { /// @ref https://en.wikipedia.org/wiki/Great-circle_navigation /// std::pair greatCircleCourse(const Point2& lonLat1, - const Point2& lonLat2) { + const Point2& lonLat2) { + if (lonLat1 == lonLat2) { + return std::make_pair(0., 0.); + } const auto lambda1 = lonLat1[0] * util::Constants::degreesToRadians(); const auto lambda2 = lonLat2[0] * util::Constants::degreesToRadians(); @@ -71,71 +70,76 @@ std::pair greatCircleCourse(const Point2& lonLat1, alpha2 * util::Constants::radiansToDegrees()); }; - -} // namespace detail -} // namespace geometry +} // namespace detail +} // namespace geometry extern "C" { // ------------------------------------------------------------------ // C wrapper interfaces to C++ routines Geometry::Implementation* atlas__Geometry__new_name(const char* name) { - Geometry::Implementation* geometry; - { - Geometry handle{std::string{name}}; - geometry = handle.get(); - geometry->attach(); - } - geometry->detach(); - return geometry; + Geometry::Implementation* geometry; + { + Geometry handle{std::string{name}}; + geometry = handle.get(); + geometry->attach(); + } + geometry->detach(); + return geometry; } -geometry::detail::GeometryBase* atlas__Geometry__new_radius(const double radius) { - Geometry::Implementation* geometry; - { - Geometry handle{radius}; - geometry = handle.get(); - geometry->attach(); - } - geometry->detach(); - return geometry; +geometry::detail::GeometryBase* atlas__Geometry__new_radius( + const double radius) { + Geometry::Implementation* geometry; + { + Geometry handle{radius}; + geometry = handle.get(); + geometry->attach(); + } + geometry->detach(); + return geometry; } void atlas__Geometry__delete(Geometry::Implementation* This) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - delete This; + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + delete This; } -void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x, const double y, const double z, - double& lon, double& lat) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - PointLonLat lonlat; - This->xyz2lonlat(PointXYZ{x, y, z}, lonlat); - lon = lonlat.lon(); - lat = lonlat.lat(); +void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x, + const double y, const double z, double& lon, + double& lat) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + PointLonLat lonlat; + This->xyz2lonlat(PointXYZ{x, y, z}, lonlat); + lon = lonlat.lon(); + lat = lonlat.lat(); } -void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This, const double lon, const double lat, double& x, +void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This, + const double lon, const double lat, double& x, double& y, double& z) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - PointXYZ xyz; - This->lonlat2xyz(PointLonLat{lon, lat}, xyz); - x = xyz.x(); - y = xyz.y(); - z = xyz.z(); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + PointXYZ xyz; + This->lonlat2xyz(PointLonLat{lon, lat}, xyz); + x = xyz.x(); + y = xyz.y(); + z = xyz.z(); } -double atlas__Geometry__distance_lonlat(Geometry::Implementation* This, const double lon1, const double lat1, +double atlas__Geometry__distance_lonlat(Geometry::Implementation* This, + const double lon1, const double lat1, const double lon2, const double lat2) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2}); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2}); } -double atlas__Geometry__distance_xyz(Geometry::Implementation* This, const double x1, const double y1, const double z1, - const double x2, const double y2, const double z2) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2}); +double atlas__Geometry__distance_xyz(Geometry::Implementation* This, + const double x1, const double y1, + const double z1, const double x2, + const double y2, const double z2) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2}); } double atlas__Geometry__radius(Geometry::Implementation* This) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->radius(); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->radius(); } double atlas__Geometry__area(Geometry::Implementation* This) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->area(); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->area(); } } // ------------------------------------------------------------------ diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index 19c778fdd..e0a1987ba 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -18,11 +18,10 @@ #include "atlas/mesh.h" #include "atlas/meshgenerator.h" #include "atlas/output/Gmsh.h" -#include "atlas/util/Point.h" #include "atlas/util/Config.h" #include "atlas/util/Constants.h" +#include "atlas/util/Point.h" #include "atlas/util/function/VortexRollup.h" - #include "tests/AtlasTestEnvironment.h" namespace atlas { @@ -37,14 +36,11 @@ constexpr auto Rank3dField = 3; // Return (u, v) field with vortex_rollup as the streamfunction. // This has no physical significance, but it makes a nice swirly field. std::pair vortexHorizontal(double lon, double lat) { - // set hLon and hLat step size. const double hLon = 0.0001; const double hLat = 0.0001; // Get finite differences. - - // Set u. const double u = (function::vortex_rollup(lon, lat + 0.5 * hLat, 0.1) - function::vortex_rollup(lon, lat - 0.5 * hLat, 0.1)) / hLat; @@ -61,7 +57,6 @@ double vortexVertical(double lon, double lat) { } void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { - const auto functionSpace = fieldSet[0].functionspace(); const auto structuredColums = functionspace::StructuredColumns(functionSpace); const auto nodeColumns = functionspace::NodeColumns(functionSpace); @@ -72,7 +67,6 @@ void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { Config("info", true); const auto gmsh = output::Gmsh(fileName, gmshConfig); - gmsh.write(mesh); gmsh.write(fieldSet, functionSpace); } @@ -118,7 +112,6 @@ struct FieldSpecFixtures { // Helper stcut to key different interpolation schemes to strings struct InterpSchemeFixtures { static const Config& get(const std::string& fixture) { - static const auto cubedsphereBilinear = option::type("cubedsphere-bilinear"); static const auto finiteElement = option::type("finite-element"); @@ -147,8 +140,9 @@ double dotProduct(const array::ArrayView& a, const array::ArrayView& b) { auto dotProd = 0.; arrayForEachDim(std::make_integer_sequence{}, std::tie(a, b), - [&](const double& aElem, - const double& bElem) { dotProd += aElem * bElem; }); + [&](const double& aElem, const double& bElem) { + dotProd += aElem * bElem; + }); return dotProd; } @@ -157,16 +151,15 @@ int countNans(const array::ArrayView& view) { auto nNans = 0; arrayForEachDim(std::make_integer_sequence{}, std::tie(view), [&](const double& viewElem) { - if (!std::isfinite(viewElem)) { - ++nNans; - } - }); + if (!std::isfinite(viewElem)) { + ++nNans; + } + }); return nNans; } template void testInterpolation(const Config& config) { - const auto& sourceFunctionSpace = FunctionSpaceFixtures::get(config.getString("source_fixture")); const auto& targetFunctionSpace = @@ -192,22 +185,22 @@ void testInterpolation(const Config& config) { auto sourceView = array::make_view(sourceField); auto targetView = array::make_view(targetField); - ArrayForEach<0>::apply(std::tie(sourceLonLat, sourceView), - [](auto&& lonLat, auto&& sourceColumn) { - - const auto setElems = [&](auto&& sourceElem) { - std::tie(sourceElem(0), sourceElem(1)) = - vortexHorizontal(lonLat(0), lonLat(1)); - if (sourceElem.size() == 3) { - sourceElem(2) = vortexVertical(lonLat(0), lonLat(1)); - } - }; - if constexpr (Rank == 2) { setElems(sourceColumn); } - else if constexpr (Rank == 3) { - ArrayForEach<0>::apply(std::tie(sourceColumn), setElems); - } - }); - sourceFieldSet.set_dirty(false); + ArrayForEach<0>::apply( + std::tie(sourceLonLat, sourceView), + [](auto&& lonLat, auto&& sourceColumn) { + const auto setElems = [&](auto&& sourceElem) { + std::tie(sourceElem(0), sourceElem(1)) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (sourceElem.size() == 3) { + sourceElem(2) = vortexVertical(lonLat(0), lonLat(1)); + } + }; + if constexpr (Rank == 2) { + setElems(sourceColumn); + } else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(sourceColumn), setElems); + } + }); const auto interp = Interpolation( InterpSchemeFixtures::get(config.getString("interp_fixture")), @@ -224,42 +217,40 @@ void testInterpolation(const Config& config) { errorView.assign(0.); auto maxError = 0.; - ArrayForEach<0>::apply(std::tie(targetLonLat, targetView, errorView), - [&](auto&& lonLat, auto&& targetColumn, - auto&& errorColumn) { - - const auto calcError = [&](auto&& targetElem, auto&& errorElem) { - auto trueValue = std::vector(targetElem.size()); - std::tie(trueValue[0], trueValue[1]) = - vortexHorizontal(lonLat(0), lonLat(1)); - if (targetElem.size() == 3) { - trueValue[2] = vortexVertical(lonLat(0), lonLat(1)); - } - - auto errorSqrd = 0.; - for (auto k = 0; k < targetElem.size(); ++k) { - errorSqrd += - (targetElem(k) - trueValue[k]) * (targetElem(k) - trueValue[k]); - } - - errorElem = std::sqrt(errorSqrd); - maxError = std::max(maxError, static_cast(errorElem)); - }; - - if - constexpr(Rank == 2) { calcError(targetColumn, errorColumn); } - else if - constexpr(Rank == 3) { - ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), calcError); - } - }); + ArrayForEach<0>::apply( + std::tie(targetLonLat, targetView, errorView), + [&](auto&& lonLat, auto&& targetColumn, auto&& errorColumn) { + const auto calcError = [&](auto&& targetElem, auto&& errorElem) { + auto trueValue = std::vector(targetElem.size()); + std::tie(trueValue[0], trueValue[1]) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (targetElem.size() == 3) { + trueValue[2] = vortexVertical(lonLat(0), lonLat(1)); + } + + auto errorSqrd = 0.; + for (auto k = 0; k < targetElem.size(); ++k) { + errorSqrd += + (targetElem(k) - trueValue[k]) * (targetElem(k) - trueValue[k]); + } + + errorElem = std::sqrt(errorSqrd); + maxError = std::max(maxError, static_cast(errorElem)); + }; + + if constexpr (Rank == 2) { + calcError(targetColumn, errorColumn); + } else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), + calcError); + } + }); EXPECT_APPROX_EQ(maxError, 0., config.getDouble("tol")); gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); - // Adjoint test auto targetAdjoint = targetFunctionSpace.createField(fieldSpec); auto targetAdjointView = array::make_view(targetAdjoint); @@ -269,8 +260,8 @@ void testInterpolation(const Config& config) { auto sourceAdjoint = sourceFunctionSpace.createField(fieldSpec); auto sourceAdjointView = array::make_view(sourceAdjoint); sourceAdjointView.assign(0.); - sourceAdjoint.set_dirty(false); + sourceAdjoint.set_dirty(false); interp.execute_adjoint(sourceAdjoint, targetAdjoint); // Check fields for nans or +/-inf @@ -311,57 +302,50 @@ CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { } CASE("finite element vector interpolation (2d-field, 2-vector)") { - const auto config = - Config("source_fixture", "gaussian_mesh") - .set("target_fixture", "cubedsphere_mesh") - .set("field_spec_fixture", "2vector") - .set("interp_fixture", "finite_element_spherical") - .set("file_id", "spherical_vector_fe") - .set("tol", 0.00015); + const auto config = Config("source_fixture", "gaussian_mesh") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "finite_element_spherical") + .set("file_id", "spherical_vector_fe") + .set("tol", 0.00015); testInterpolation((config)); } CASE("structured columns vector interpolation (2d-field, 2-vector)") { - - const auto config = - Config("source_fixture", "structured_columns") - .set("target_fixture", "cubedsphere_mesh") - .set("field_spec_fixture", "2vector") - .set("interp_fixture", "structured_linear_spherical") - .set("file_id", "spherical_vector_sc") - .set("tol", 0.00017); + const auto config = Config("source_fixture", "structured_columns") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc") + .set("tol", 0.00017); testInterpolation((config)); } CASE("structured columns vector interpolation (2d-field, 2-vector, low-res)") { - - const auto config = - Config("source_fixture", "structured_columns_lowres") - .set("target_fixture", "gaussian_mesh") - .set("field_spec_fixture", "2vector") - .set("interp_fixture", "structured_linear_spherical") - .set("file_id", "spherical_vector_sc_lr") - .set("tol", 0.00056); + const auto config = Config("source_fixture", "structured_columns_lowres") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_lr") + .set("tol", 0.00056); testInterpolation((config)); } CASE("structured columns vector interpolation (2d-field, 2-vector, hi-res)") { - - const auto config = - Config("source_fixture", "structured_columns_hires") - .set("target_fixture", "gaussian_mesh") - .set("field_spec_fixture", "2vector") - .set("interp_fixture", "structured_linear_spherical") - .set("file_id", "spherical_vector_sc_hr") - .set("tol", 0.000044); + const auto config = Config("source_fixture", "structured_columns_hires") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_hr") + .set("tol", 0.000044); testInterpolation((config)); } -} -} +} // namespace test +} // namespace atlas int main(int argc, char** argv) { return atlas::test::run(argc, argv); } diff --git a/src/tests/util/test_unitsphere.cc b/src/tests/util/test_unitsphere.cc index d390334a5..8bc249c6a 100644 --- a/src/tests/util/test_unitsphere.cc +++ b/src/tests/util/test_unitsphere.cc @@ -5,9 +5,8 @@ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ -#include "atlas/util/Point.h" #include "atlas/util/Geometry.h" - +#include "atlas/util/Point.h" #include "tests/AtlasTestEnvironment.h" namespace atlas { @@ -16,28 +15,43 @@ namespace test { using namespace atlas::util; CASE("great-circle course") { - geometry::UnitSphere g; const auto point1 = PointLonLat(-71.6, -33.0); // Valparaiso const auto point2 = PointLonLat(121.8, 31.4); // Shanghai const auto point3 = PointLonLat(0., 89.); const auto point4 = PointLonLat(180., 89.); + const auto point5 = PointLonLat(0., 90.); + const auto point6 = PointLonLat(180., 90.); const auto targetCourse1 = -94.41; const auto targetCourse2 = -78.42; const auto targetCourse3 = 0.; const auto targetCourse4 = 180.; + const auto targetCourse5 = 0.; + const auto targetCourse6 = 180.; + const auto targetCourse7 = 0.; + const auto targetCourse8 = 0.; + + + const auto [course1, course2] = g.greatCircleCourse(point1, point2); + const auto [course3, course4] = g.greatCircleCourse(point3, point4); + + // Colocated points on pole. + const auto [course5, course6] = g.greatCircleCourse(point5, point6); + const auto [course7, course8] = g.greatCircleCourse(point5, point5); - const auto[ course1, course2 ] = g.greatCircleCourse(point1, point2); - const auto[ course3, course4 ] = g.greatCircleCourse(point3, point4); EXPECT_APPROX_EQ(course1, targetCourse1, 0.01); EXPECT_APPROX_EQ(course2, targetCourse2, 0.01); EXPECT_APPROX_EQ(course3, targetCourse3, 1.e-14); EXPECT_APPROX_EQ(std::abs(course4), targetCourse4, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course5), targetCourse5, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course6), targetCourse6, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course7), targetCourse7, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course8), targetCourse8, 1.e-14); } -} // namespace test -} // namespace atlas +} // namespace test +} // namespace atlas int main(int argc, char** argv) { return atlas::test::run(argc, argv); } From de3896fa4d3102fbe2a69819174978eba7b11087 Mon Sep 17 00:00:00 2001 From: Oliver Lomax Date: Wed, 6 Mar 2024 22:27:00 +0000 Subject: [PATCH 12/23] Removed leftover code missed in PR #175 --- .../method/sphericalvector/SphericalVector.cc | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index a8cc3a96f..7e29f41c1 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -72,15 +72,9 @@ void SphericalVector::do_setup(const FunctionSpace& source, // whereas eckit does not. auto complexTriplets = ComplexTriplets(nNonZeros); auto realTriplets = RealTriplets(nNonZeros); - - // Make sure halo lonlats are same as owned points. - auto sourceLonLats = source_.createField(option::name("lonlat") | - option::variables(2)); - auto targetLonLats = target_.createField(option::name("lonlat") | - option::variables(2)); + const auto sourceLonLatsView = array::make_view(source_.lonlat()); const auto targetLonLatsView = array::make_view(target_.lonlat()); - const auto unitSphere = geometry::UnitSphere{}; atlas_omp_parallel_for(auto rowIndex = Index{0}; rowIndex < nRows; From 5dc85c9afe3f256787ed775f19288671db03c2c0 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 25 Jan 2024 09:04:29 +0100 Subject: [PATCH 13/23] Use new LocalConfiguration baseclass functions in util::Config and util::Metadata instead of eckit::Value backdoor This is available since ecmwf/eckit#105 --- atlas_io/CMakeLists.txt | 4 ++++ atlas_io/src/atlas_io/CMakeLists.txt | 11 +++++++++++ atlas_io/src/atlas_io/Metadata.h | 16 ++++++++++++++-- atlas_io/src/atlas_io/detail/defines.h.in | 5 +++++ src/atlas/util/Config.cc | 18 ++++++++++++++---- src/atlas/util/Config.h | 8 ++++++-- src/atlas/util/Metadata.cc | 6 +++++- src/atlas/util/Metadata.h | 2 +- 8 files changed, 60 insertions(+), 10 deletions(-) diff --git a/atlas_io/CMakeLists.txt b/atlas_io/CMakeLists.txt index 26f549808..8417ee369 100644 --- a/atlas_io/CMakeLists.txt +++ b/atlas_io/CMakeLists.txt @@ -28,6 +28,10 @@ ecbuild_debug( " eckit_FEATURES : [${eckit_FEATURES}]" ) ################################################################################ # Features that can be enabled / disabled with -DENABLE_ +ecbuild_add_option( FEATURE ECKIT_DEVELOP + DESCRIPTION "Used to enable new features or API depending on eckit develop branch, not yet in a tagged release" + DEFAULT OFF ) + ################################################################################ # sources diff --git a/atlas_io/src/atlas_io/CMakeLists.txt b/atlas_io/src/atlas_io/CMakeLists.txt index 877d61ead..6113287b2 100644 --- a/atlas_io/src/atlas_io/CMakeLists.txt +++ b/atlas_io/src/atlas_io/CMakeLists.txt @@ -1,3 +1,14 @@ + + +if( HAVE_ECKIT_DEVELOP ) + set( ATLAS_IO_ECKIT_DEVELOP 1 ) +else() + set( ATLAS_IO_ECKIT_DEVELOP 0 ) +endif() + +ecbuild_parse_version( ${eckit_VERSION} PREFIX ATLAS_IO_ECKIT ) +math( EXPR ATLAS_IO_ECKIT_VERSION_INT "( 10000 * ${ATLAS_IO_ECKIT_VERSION_MAJOR} ) + ( 100 * ${ATLAS_IO_ECKIT_VERSION_MINOR} ) + ${ATLAS_IO_ECKIT_VERSION_PATCH}" ) + configure_file( detail/defines.h.in detail/defines.h ) install( FILES ${CMAKE_CURRENT_BINARY_DIR}/detail/defines.h diff --git a/atlas_io/src/atlas_io/Metadata.h b/atlas_io/src/atlas_io/Metadata.h index b6f2ab123..d6689ba21 100644 --- a/atlas_io/src/atlas_io/Metadata.h +++ b/atlas_io/src/atlas_io/Metadata.h @@ -57,7 +57,11 @@ class Metadata : public eckit::LocalConfiguration { // extended LocalConfiguration: using eckit::LocalConfiguration::set; - Metadata& set(const eckit::LocalConfiguration& other) { + + Metadata& set(const eckit::Configuration& other) { +#if ATLAS_IO_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_IO_ECKIT_DEVELOP + LocalConfiguration::set(other); +#else eckit::Value& root = const_cast(get()); auto& other_root = other.get(); std::vector other_keys; @@ -65,6 +69,7 @@ class Metadata : public eckit::LocalConfiguration { for (auto& key : other_keys) { root[key] = other_root[key]; } +#endif return *this; } @@ -79,17 +84,24 @@ class Metadata : public eckit::LocalConfiguration { Metadata& remove(const std::string& name) { +#if ATLAS_IO_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_IO_ECKIT_DEVELOP + LocalConfiguration::remove(name); +#else eckit::Value& root = const_cast(get()); root.remove(name); +#endif return *this; } std::vector keys() const { - // Preserves order of keys +#if ATLAS_IO_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_IO_ECKIT_DEVELOP + return LocalConfiguration::keys(); +#else std::vector result; eckit::fromValue(result, get().keys()); return result; +#endif } }; diff --git a/atlas_io/src/atlas_io/detail/defines.h.in b/atlas_io/src/atlas_io/detail/defines.h.in index 40d8132c8..f3fb79f42 100644 --- a/atlas_io/src/atlas_io/detail/defines.h.in +++ b/atlas_io/src/atlas_io/detail/defines.h.in @@ -18,4 +18,9 @@ #cmakedefine01 ATLAS_IO_LITTLE_ENDIAN #cmakedefine01 ATLAS_IO_BIG_ENDIAN +#define ATLAS_IO_ECKIT_VERSION_INT @ATLAS_IO_ECKIT_VERSION_INT@ +#define ATLAS_IO_ECKIT_DEVELOP @ATLAS_IO_ECKIT_DEVELOP@ + +#define ATLAS_IO_ECKIT_VERSION_AT_LEAST(x, y, z) (ATLAS_IO_ECKIT_VERSION_INT >= x * 10000 + y * 100 + z) + #endif diff --git a/src/atlas/util/Config.cc b/src/atlas/util/Config.cc index 8ccf234c1..b54ba2318 100644 --- a/src/atlas/util/Config.cc +++ b/src/atlas/util/Config.cc @@ -59,13 +59,16 @@ Config::Config(std::istream& stream, const std::string&): eckit::LocalConfigurat Config::Config(const eckit::PathName& path): eckit::LocalConfiguration(yaml_from_path(path)) {} -Config Config::operator|(const Config& other) const { +Config Config::operator|(const eckit::Configuration& other) const { Config config(*this); config.set(other); return config; } -Config& Config::set(const eckit::LocalConfiguration& other) { +Config& Config::set(const eckit::Configuration& other) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP + eckit::LocalConfiguration::set(other); +#else eckit::Value& root = const_cast(get()); auto& other_root = other.get(); std::vector other_keys; @@ -73,12 +76,17 @@ Config& Config::set(const eckit::LocalConfiguration& other) { for (auto& key : other_keys) { root[key] = other_root[key]; } +#endif return *this; } Config& Config::remove(const std::string& name) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP + LocalConfiguration::remove(name); +#else eckit::Value& root = const_cast(get()); root.remove(name); +#endif return *this; } @@ -94,20 +102,22 @@ Config& Config::set(const std::string& name, const std::vector& values) bool Config::get(const std::string& name, std::vector& value) const { bool found = has(name); if (found) { - std::vector properties = getSubConfigurations(name); + auto properties = getSubConfigurations(name); value.resize(properties.size()); for (size_t i = 0; i < value.size(); ++i) { - value[i] = Config(properties[i]); + value[i].set(properties[i]); } } return found; } +#if ! (ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP) std::vector Config::keys() const { std::vector result; eckit::fromValue(result, get().keys()); return result; } +#endif std::string Config::json(eckit::JSON::Formatting formatting) const { std::stringstream json; diff --git a/src/atlas/util/Config.h b/src/atlas/util/Config.h index 6c9e9a42b..5032a17f2 100644 --- a/src/atlas/util/Config.h +++ b/src/atlas/util/Config.h @@ -16,6 +16,8 @@ #include "eckit/config/LocalConfiguration.h" #include "eckit/log/JSON.h" +#include "atlas/library/config.h" + namespace eckit { class PathName; class Hash; @@ -61,14 +63,14 @@ class Config : public eckit::LocalConfiguration { Config operator()(const std::string& name, std::initializer_list&& value); // Overload operators to merge two Config objects. - Config operator|(const Config& other) const; + Config operator|(const eckit::Configuration& other) const; /// @brief Set a key-value parameter using eckit::LocalConfiguration::set; Config& set(const std::string& name, const std::vector&); - Config& set(const eckit::LocalConfiguration&); + Config& set(const eckit::Configuration&); template Config& set(const std::string& name, std::initializer_list&& value); @@ -80,7 +82,9 @@ class Config : public eckit::LocalConfiguration { using eckit::LocalConfiguration::get; bool get(const std::string& name, std::vector& value) const; +#if ! (ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP) std::vector keys() const; +#endif std::string json(eckit::JSON::Formatting = eckit::JSON::Formatting::indent()) const; }; diff --git a/src/atlas/util/Metadata.cc b/src/atlas/util/Metadata.cc index c24e6c897..583332bec 100644 --- a/src/atlas/util/Metadata.cc +++ b/src/atlas/util/Metadata.cc @@ -119,7 +119,10 @@ void Metadata::broadcast(Metadata& dest, idx_t root) const { } -Metadata& Metadata::set(const eckit::LocalConfiguration& other) { +Metadata& Metadata::set(const eckit::Configuration& other) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP + LocalConfiguration::set(other); +#else eckit::Value& root = const_cast(get()); auto& other_root = other.get(); std::vector other_keys; @@ -127,6 +130,7 @@ Metadata& Metadata::set(const eckit::LocalConfiguration& other) { for (auto& key : other_keys) { root[key] = other_root[key]; } +#endif return *this; } diff --git a/src/atlas/util/Metadata.h b/src/atlas/util/Metadata.h index c697e13fa..f4ce34eef 100644 --- a/src/atlas/util/Metadata.h +++ b/src/atlas/util/Metadata.h @@ -35,7 +35,7 @@ class Metadata : public eckit::LocalConfiguration { return *this; } - Metadata& set(const eckit::LocalConfiguration&); + Metadata& set(const eckit::Configuration&); using eckit::LocalConfiguration::get; From 8b37fca62913c6429471df27f5d84a48d45c3f3a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 21 Mar 2024 16:47:46 +0100 Subject: [PATCH 14/23] Avoid linker warnings on macOS about 'ld: warning: could not create compact unwind for ...' --- cmake/atlas_compile_flags.cmake | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cmake/atlas_compile_flags.cmake b/cmake/atlas_compile_flags.cmake index 8b334f149..5a2274f08 100644 --- a/cmake/atlas_compile_flags.cmake +++ b/cmake/atlas_compile_flags.cmake @@ -32,3 +32,7 @@ if( CMAKE_CXX_COMPILER_ID MATCHES IntelLLVM ) ecbuild_add_cxx_flags("-fno-finite-math-only") endif() +if( APPLE AND CMAKE_Fortran_COMPILER_ID MATCHES GNU ) + # Avoid warnings : ld: warning: could not create compact unwind for ... + add_link_options(-Wl,-keep_dwarf_unwind -Wl,-no_compact_unwind) +endif() From 56a7599ed4866d95f61e5c0985e0ded17e3d9655 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 22 Mar 2024 10:59:22 +0000 Subject: [PATCH 15/23] Fix bug in TraceT caused by typo where the title was wrong --- src/atlas/runtime/trace/TraceT.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/atlas/runtime/trace/TraceT.h b/src/atlas/runtime/trace/TraceT.h index efe98f45c..21def5df8 100644 --- a/src/atlas/runtime/trace/TraceT.h +++ b/src/atlas/runtime/trace/TraceT.h @@ -93,8 +93,8 @@ class TraceT { template inline std::string TraceT::formatTitle(const std::string& _title) { - std::string title = _title; - +(Barriers::state() ? " [b]" : "") + + std::string title = + _title + (Barriers::state() ? " [b]" : "") + (atlas_omp_get_num_threads() > 1 ? " @thread[" + std::to_string(atlas_omp_get_thread_num()) + "]" : ""); return title; } From a7d336caee35e389ac7a0514d511d39db702ee25 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 22 Mar 2024 10:59:36 +0000 Subject: [PATCH 16/23] Initialize std::array values to zero because valgrind complains, even though c++ standard mandates it should be default-initialized to zero --- src/atlas/output/detail/GmshIO.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/atlas/output/detail/GmshIO.cc b/src/atlas/output/detail/GmshIO.cc index 305051de4..810c30ec1 100644 --- a/src/atlas/output/detail/GmshIO.cc +++ b/src/atlas/output/detail/GmshIO.cc @@ -143,6 +143,7 @@ void write_level(std::ostream& out, GlobalIndex gidx, const array::LocalView data_vec; + data_vec.fill(static_cast(0)); for (idx_t n = 0; n < ndata; ++n) { if (include(n)) { for (idx_t v = 0; v < nvars; ++v) { @@ -158,6 +159,7 @@ void write_level(std::ostream& out, GlobalIndex gidx, const array::LocalView data_vec; + data_vec.fill(static_cast(0)); if (nvars == 4) { for (int n = 0; n < ndata; ++n) { if (include(n)) { From 0f785d1fd2d089844b8e6958e10ef2bd3082d319 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 22 Mar 2024 10:59:58 +0000 Subject: [PATCH 17/23] Cosmetic: readability braces --- .../interpolation/test_interpolation_spherical_vector.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index e0a1987ba..4717c5a19 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -175,7 +175,9 @@ void testInterpolation(const Config& config) { auto fieldSpec = FieldSpecFixtures::get(config.getString("field_spec_fixture")); - if constexpr (Rank == 3) fieldSpec.set("levels", 2); + if constexpr (Rank == 3) { + fieldSpec.set("levels", 2); + } auto sourceField = sourceFieldSet.add(sourceFunctionSpace.createField(fieldSpec)); From 7e1134a3752f8d5b57bf318363a536efb247e787 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 22 Mar 2024 14:07:04 +0100 Subject: [PATCH 18/23] Revert "Avoid linker warnings on macOS about 'ld: warning: could not create compact unwind for ...'" This reverts commit 8b37fca62913c6429471df27f5d84a48d45c3f3a. --- cmake/atlas_compile_flags.cmake | 4 ---- 1 file changed, 4 deletions(-) diff --git a/cmake/atlas_compile_flags.cmake b/cmake/atlas_compile_flags.cmake index 5a2274f08..8b334f149 100644 --- a/cmake/atlas_compile_flags.cmake +++ b/cmake/atlas_compile_flags.cmake @@ -32,7 +32,3 @@ if( CMAKE_CXX_COMPILER_ID MATCHES IntelLLVM ) ecbuild_add_cxx_flags("-fno-finite-math-only") endif() -if( APPLE AND CMAKE_Fortran_COMPILER_ID MATCHES GNU ) - # Avoid warnings : ld: warning: could not create compact unwind for ... - add_link_options(-Wl,-keep_dwarf_unwind -Wl,-no_compact_unwind) -endif() From 3bcbdb8ef0bf4d0f6befc653ded3cb37851829f8 Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Mon, 25 Mar 2024 10:45:25 +0000 Subject: [PATCH 19/23] Fix build for configuration setting ATLAS_BITS_LOCAL=64 (#184) --- 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 61388a2f5..702e4bef5 100644 --- a/src/atlas/functionspace/PointCloud.cc +++ b/src/atlas/functionspace/PointCloud.cc @@ -1001,7 +1001,7 @@ void PointCloud::setupGatherScatter() { array::make_view(remote_index_).data(), REMOTE_IDX_BASE, array::make_view(global_index_).data(), - array::make_view(ghost_).data(), + array::make_view(ghost_).data(), ghost_.size()); size_global_ = gather_scatter_->glb_dof(); } From 7abbf1a41e1aa207a94243435f8493a9d14562e4 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 25 Mar 2024 17:13:46 +0100 Subject: [PATCH 20/23] atlas_io is an adaptor library when eckit_codec is available (#181) --- atlas_io/CMakeLists.txt | 18 +++- atlas_io/eckit_codec_adaptor/CMakeLists.txt | 1 + .../eckit_codec_adaptor/src/CMakeLists.txt | 2 + .../src/atlas_io/CMakeLists.txt | 19 ++++ .../eckit_codec_adaptor/src/atlas_io/Trace.cc | 43 ++++++++ .../eckit_codec_adaptor/src/atlas_io/Trace.h | 77 +++++++++++++ .../src/atlas_io/atlas-io.h | 99 +++++++++++++++++ .../src/atlas_io/detail/BlackMagic.h | 101 ++++++++++++++++++ atlas_io/tests/TestEnvironment.h | 2 +- atlas_io/tests/test_io_stream.cc | 4 +- 10 files changed, 360 insertions(+), 6 deletions(-) create mode 100644 atlas_io/eckit_codec_adaptor/CMakeLists.txt create mode 100644 atlas_io/eckit_codec_adaptor/src/CMakeLists.txt create mode 100644 atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt create mode 100644 atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc create mode 100644 atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h create mode 100644 atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h create mode 100644 atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h diff --git a/atlas_io/CMakeLists.txt b/atlas_io/CMakeLists.txt index 8417ee369..2a821b761 100644 --- a/atlas_io/CMakeLists.txt +++ b/atlas_io/CMakeLists.txt @@ -32,6 +32,16 @@ ecbuild_add_option( FEATURE ECKIT_DEVELOP DESCRIPTION "Used to enable new features or API depending on eckit develop branch, not yet in a tagged release" DEFAULT OFF ) +set( eckit_HAVE_ECKIT_CODEC 0 ) +if( TARGET eckit_codec ) + set( eckit_HAVE_ECKIT_CODEC 1 ) +endif() + +ecbuild_add_option( FEATURE ECKIT_CODEC + DEFAULT ON + DESCRIPTION "Use eckit::codec with adaptor" + CONDITION eckit_HAVE_ECKIT_CODEC ) + ################################################################################ # sources @@ -51,8 +61,12 @@ else() set( ATLAS_IO_LITTLE_ENDIAN 1 ) endif() - -add_subdirectory( src ) +if( HAVE_ECKIT_CODEC ) + ecbuild_info("atlas_io is configured to be an adaptor library which delegates calls to eckit_codec") + add_subdirectory(eckit_codec_adaptor) +else() + add_subdirectory( src ) +endif() add_subdirectory( tests ) ################################################################################ diff --git a/atlas_io/eckit_codec_adaptor/CMakeLists.txt b/atlas_io/eckit_codec_adaptor/CMakeLists.txt new file mode 100644 index 000000000..febd4f0ab --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(src) diff --git a/atlas_io/eckit_codec_adaptor/src/CMakeLists.txt b/atlas_io/eckit_codec_adaptor/src/CMakeLists.txt new file mode 100644 index 000000000..e91bc52fd --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(atlas_io) + diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt b/atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt new file mode 100644 index 000000000..1794974b1 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt @@ -0,0 +1,19 @@ + +ecbuild_add_library( TARGET atlas_io + + INSTALL_HEADERS ALL + HEADER_DESTINATION include/atlas_io + PUBLIC_LIBS eckit_codec + PUBLIC_INCLUDES + $ + + SOURCES + atlas-io.h + Trace.cc + Trace.h + detail/BlackMagic.h + LINKER_LANGUAGE CXX +) + +target_compile_features( atlas_io PUBLIC cxx_std_17 ) + diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc new file mode 100644 index 000000000..5e52b7dc9 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc @@ -0,0 +1,43 @@ +/* + * (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. + */ + +#include "Trace.h" + +#include "eckit/log/CodeLocation.h" + +namespace atlas { +namespace io { + +atlas::io::Trace::Trace(const eckit::CodeLocation& loc) { + for (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, loc.func())); + } + } +} + +Trace::Trace(const eckit::CodeLocation& loc, const std::string& title) { + for (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, title)); + } + } +} + +Trace::Trace(const eckit::CodeLocation& loc, const std::string& title, const Labels& labels) { + for (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, title)); + } + } +} + +} // namespace io +} // namespace atlas diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h new file mode 100644 index 000000000..43da3f526 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h @@ -0,0 +1,77 @@ +/* + * (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 +#include +#include + +namespace eckit { +class CodeLocation; +} + +namespace atlas { +namespace io { + +struct TraceHook { + TraceHook() = default; + virtual ~TraceHook() = default; +}; + +struct TraceHookRegistry { + using TraceHookBuilder = std::function(const eckit::CodeLocation&, const std::string&)>; + std::vector hooks; + std::vector enabled_; + static TraceHookRegistry& instance() { + static TraceHookRegistry instance; + return instance; + } + static size_t add(TraceHookBuilder&& hook) { + instance().hooks.emplace_back(hook); + instance().enabled_.emplace_back(true); + return instance().hooks.size() - 1; + } + static size_t add(const TraceHookBuilder& hook) { + instance().hooks.emplace_back(hook); + instance().enabled_.emplace_back(true); + return instance().hooks.size() - 1; + } + static void enable(size_t id) { instance().enabled_[id] = true; } + static void disable(size_t id) { instance().enabled_[id] = false; } + static bool enabled(size_t id) { return instance().enabled_[id]; } + static size_t size() { return instance().hooks.size(); } + static TraceHookBuilder& hook(size_t id) { return instance().hooks[id]; } + static size_t invalidId() { return std::numeric_limits::max(); } + +private: + TraceHookRegistry() = default; +}; + +struct Trace { + using Labels = std::vector; + Trace(const eckit::CodeLocation& loc); + Trace(const eckit::CodeLocation& loc, const std::string& title); + Trace(const eckit::CodeLocation& loc, const std::string& title, const Labels& labels); + +private: + std::vector> hooks_; +}; + +} // namespace io +} // namespace atlas + +#include "atlas_io/detail/BlackMagic.h" + +#define ATLAS_IO_TRACE(...) __ATLAS_IO_TYPE(::atlas::io::Trace, Here() __ATLAS_IO_COMMA_ARGS(__VA_ARGS__)) +#define ATLAS_IO_TRACE_SCOPE(...) __ATLAS_IO_TYPE_SCOPE(::atlas::io::Trace, Here() __ATLAS_IO_COMMA_ARGS(__VA_ARGS__)) diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h b/atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h new file mode 100644 index 000000000..7c5f49d20 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h @@ -0,0 +1,99 @@ +/* + * (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 "eckit/codec/codec.h" + +namespace atlas::io { + + // Encoding/Decoding + using Metadata = ::eckit::codec::Metadata; + using DataType = ::eckit::codec::DataType; + using Data = ::eckit::codec::Data; + using Encoder = ::eckit::codec::Encoder; + using Decoder = ::eckit::codec::Decoder; + + // Record + using Record = ::eckit::codec::Record; + using RecordWriter = ::eckit::codec::RecordWriter; + using RecordPrinter = ::eckit::codec::RecordPrinter; + using RecordItemReader = ::eckit::codec::RecordItemReader; + using RecordReader = ::eckit::codec::RecordReader; + + // I/O + using Session = ::eckit::codec::Session; + using Mode = ::eckit::codec::Mode; + using Stream = ::eckit::codec::Stream; + using FileStream = ::eckit::codec::FileStream; + using InputFileStream = ::eckit::codec::InputFileStream; + using OutputFileStream = ::eckit::codec::OutputFileStream; + + // Array + using ArrayReference = ::eckit::codec::ArrayReference; + using ArrayMetadata = ::eckit::codec::ArrayMetadata; + using ArrayShape = ::eckit::codec::ArrayShape; + + // Exceptions + using Exception = ::eckit::codec::Exception; + using NotEncodable = ::eckit::codec::NotEncodable; + using NotDecodable = ::eckit::codec::NotDecodable; + using InvalidRecord = ::eckit::codec::InvalidRecord; + using DataCorruption = ::eckit::codec::DataCorruption; + using WriteError = ::eckit::codec::WriteError; + + // Type traits + template + static constexpr bool is_interpretable() { + return ::eckit::codec::is_interpretable(); + } + template + static constexpr bool is_encodable() { + return ::eckit::codec::is_encodable(); + } + template + static constexpr bool is_decodable() { + return ::eckit::codec::is_decodable(); + } + template + static constexpr bool can_encode_metadata() { + return ::eckit::codec::can_encode_metadata(); + } + template + static constexpr bool can_encode_data() { + return ::eckit::codec::can_encode_metadata(); + } + + namespace tag { + using disable_static_assert = ::eckit::codec::tag::disable_static_assert; + using enable_static_assert = ::eckit::codec::tag::enable_static_assert; + } + + // Functions + using ::eckit::codec::ref; + using ::eckit::codec::copy; + using ::eckit::codec::encode; + using ::eckit::codec::decode; + using ::eckit::codec::interprete; + using ::eckit::codec::link; + using ::eckit::codec::make_datatype; + // template + // using make_datatype = eckit::codec::make_datatype; + + namespace defaults { + using ::eckit::codec::defaults::compression_algorithm; + using ::eckit::codec::defaults::checksum_algorithm; + using ::eckit::codec::defaults::checksum_read; + using ::eckit::codec::defaults::checksum_write; + } +} + +#define ATLAS_IO_ASSERT(X) ASSERT(X) +#define ATLAS_IO_ASSERT_MSG(X, M) ASSERT_MSG(X, M) diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h b/atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h new file mode 100644 index 000000000..a889c5e19 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h @@ -0,0 +1,101 @@ +/* + * (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. + */ + +#pragma once + +// This file contains preprocessor black magic. It contains macros that +// can return the number of arguments passed + +//----------------------------------------------------------------------------------------------------------- +// Public + +/// Returns the number of passed arguments +#define __ATLAS_IO_NARG(...) + +/// Splice a and b together +#define __ATLAS_IO_SPLICE(a, b) + +#define __ATLAS_IO_STRINGIFY(a) a + +#define __ATLAS_IO_TYPE(Type, ...) +#define __ATLAS_IO_TYPE_SCOPE(Type, ...) + +//----------------------------------------------------------------------------------------------------------- +// Details + +// Undefine these, to be redefined further down. +#undef __ATLAS_IO_NARG +#undef __ATLAS_IO_SPLICE +#undef __ATLAS_IO_TYPE +#undef __ATLAS_IO_TYPE_SCOPE + +#define __ATLAS_IO_REVERSE 5, 4, 3, 2, 1, 0 +#define __ATLAS_IO_ARGN(_1, _2, _3, _4, _5, N, ...) N +#define __ATLAS_IO_NARG__(dummy, ...) __ATLAS_IO_ARGN(__VA_ARGS__) +#define __ATLAS_IO_NARG_(...) __ATLAS_IO_NARG__(dummy, ##__VA_ARGS__, __ATLAS_IO_REVERSE) +#define __ATLAS_IO_SPLICE(a, b) __ATLAS_IO_SPLICE_1(a, b) +#define __ATLAS_IO_SPLICE_1(a, b) __ATLAS_IO_SPLICE_2(a, b) +#define __ATLAS_IO_SPLICE_2(a, b) a##b + +#define __ATLAS_IO_ARG16(_0, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, _15, ...) _15 +#define __ATLAS_IO_HAS_COMMA(...) __ATLAS_IO_ARG16(__VA_ARGS__, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0) +#define __ATLAS_IO_TRIGGER_PARENTHESIS(...) , +#define __ATLAS_IO_ISEMPTY(...) \ + __ATLAS_IO_ISEMPTY_(/* test if there is just one argument, eventually an empty \ + one */ \ + __ATLAS_IO_HAS_COMMA(__VA_ARGS__), /* test if \ + _TRIGGER_PARENTHESIS_ \ + together with the \ + argument adds a comma */ \ + __ATLAS_IO_HAS_COMMA( \ + __ATLAS_IO_TRIGGER_PARENTHESIS __VA_ARGS__), /* test if the argument together with \ + a parenthesis adds a comma */ \ + __ATLAS_IO_HAS_COMMA(__VA_ARGS__(/*empty*/)), /* test if placing it between \ + __ATLAS_IO_TRIGGER_PARENTHESIS and the \ + parenthesis adds a comma */ \ + __ATLAS_IO_HAS_COMMA(__ATLAS_IO_TRIGGER_PARENTHESIS __VA_ARGS__(/*empty*/))) + +#define __ATLAS_IO_PASTE5(_0, _1, _2, _3, _4) _0##_1##_2##_3##_4 +#define __ATLAS_IO_ISEMPTY_(_0, _1, _2, _3) \ + __ATLAS_IO_HAS_COMMA(__ATLAS_IO_PASTE5(__ATLAS_IO_IS_EMPTY_CASE_, _0, _1, _2, _3)) +#define __ATLAS_IO_IS_EMPTY_CASE_0001 , + +#define __ATLAS_IO_NARG(...) __ATLAS_IO_SPLICE(__ATLAS_IO_CALL_NARG_, __ATLAS_IO_ISEMPTY(__VA_ARGS__))(__VA_ARGS__) +#define __ATLAS_IO_CALL_NARG_1(...) 0 +#define __ATLAS_IO_CALL_NARG_0 __ATLAS_IO_NARG_ + +#define __ATLAS_IO_COMMA_ARGS(...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_COMMA_ARGS_, __ATLAS_IO_ISEMPTY(__VA_ARGS__))(__VA_ARGS__) +#define __ATLAS_IO_COMMA_ARGS_1(...) +#define __ATLAS_IO_COMMA_ARGS_0(...) , __VA_ARGS__ + +#define __ATLAS_IO_ARGS_OR_DUMMY(...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_ARGS_OR_DUMMY_, __ATLAS_IO_ISEMPTY(__VA_ARGS__)) \ + (__VA_ARGS__) +#define __ATLAS_IO_ARGS_OR_DUMMY_0(...) __VA_ARGS__ +#define __ATLAS_IO_ARGS_OR_DUMMY_1(...) 0 + +#define __ATLAS_IO_TYPE(Type, ...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_TYPE_, __ATLAS_IO_ISEMPTY(__VA_ARGS__)) \ + (Type, __ATLAS_IO_ARGS_OR_DUMMY(__VA_ARGS__)) +#define __ATLAS_IO_TYPE_1(Type, dummy) Type __ATLAS_IO_SPLICE(__variable_, __LINE__) +#define __ATLAS_IO_TYPE_0(Type, ...) Type __ATLAS_IO_SPLICE(__variable_, __LINE__)(__VA_ARGS__) + +#define __ATLAS_IO_TYPE_SCOPE(Type, ...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_TYPE_SCOPE_, __ATLAS_IO_ISEMPTY(__VA_ARGS__)) \ + (Type, __ATLAS_IO_ARGS_OR_DUMMY(__VA_ARGS__)) +#define __ATLAS_IO_TYPE_SCOPE_1(Type, ...) \ + for (bool __ATLAS_IO_SPLICE(__done_, __LINE__) = false; __ATLAS_IO_SPLICE(__done_, __LINE__) != true;) \ + for (Type __ATLAS_IO_SPLICE(__variable_, __LINE__); __ATLAS_IO_SPLICE(__done_, __LINE__) != true; \ + __ATLAS_IO_SPLICE(__done_, __LINE__) = true) +#define __ATLAS_IO_TYPE_SCOPE_0(Type, ...) \ + for (bool __ATLAS_IO_SPLICE(__done_, __LINE__) = false; __ATLAS_IO_SPLICE(__done_, __LINE__) != true;) \ + for (Type __ATLAS_IO_SPLICE(__variable_, __LINE__)(__VA_ARGS__); __ATLAS_IO_SPLICE(__done_, __LINE__) != true; \ + __ATLAS_IO_SPLICE(__done_, __LINE__) = true) diff --git a/atlas_io/tests/TestEnvironment.h b/atlas_io/tests/TestEnvironment.h index aa2b99e90..7b66cee8f 100644 --- a/atlas_io/tests/TestEnvironment.h +++ b/atlas_io/tests/TestEnvironment.h @@ -27,7 +27,7 @@ #include "eckit/types/Types.h" #include "atlas_io/atlas-io.h" -#include "atlas_io/detail/BlackMagic.h" +#include "atlas_io/Trace.h" namespace atlas { namespace test { diff --git a/atlas_io/tests/test_io_stream.cc b/atlas_io/tests/test_io_stream.cc index 46ed49874..598e26b04 100644 --- a/atlas_io/tests/test_io_stream.cc +++ b/atlas_io/tests/test_io_stream.cc @@ -8,12 +8,10 @@ * nor does it submit to any jurisdiction. */ -#include "atlas_io/FileStream.h" -#include "atlas_io/Session.h" +#include "atlas_io/atlas-io.h" #include "eckit/io/FileHandle.h" #include "eckit/io/PooledHandle.h" - #include "TestEnvironment.h" From 270b936d94d111f48b2b93ff4c01bee7c3b661cf Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Tue, 9 Apr 2024 14:45:24 +0000 Subject: [PATCH 21/23] =?UTF-8?q?Projection=20base=20implementation=20deri?= =?UTF-8?q?vatives=20performance/encapsulation=20=E2=80=A6=20(#185)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Projection base implementation derivatives performance/encapsulation improvements (better use of temporaries) * Projection base implementation speedup on calculation of bounding boxes; MercatorProjection fix on handling extrema * MercatorProjection simpler handling of extrema --- src/atlas/mesh/Mesh.cc | 24 +++--- .../projection/detail/MercatorProjection.cc | 17 ++++- src/atlas/projection/detail/ProjectionImpl.cc | 76 ++++++++----------- src/atlas/projection/detail/ProjectionImpl.h | 21 ++--- 4 files changed, 63 insertions(+), 75 deletions(-) diff --git a/src/atlas/mesh/Mesh.cc b/src/atlas/mesh/Mesh.cc index 411260617..1526a74c4 100644 --- a/src/atlas/mesh/Mesh.cc +++ b/src/atlas/mesh/Mesh.cc @@ -9,11 +9,13 @@ */ #include "atlas/mesh/Mesh.h" -#include "atlas/mesh/Nodes.h" + #include "atlas/grid/Grid.h" #include "atlas/grid/Partitioner.h" +#include "atlas/mesh/Nodes.h" #include "atlas/meshgenerator/MeshGenerator.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/runtime/Exception.h" namespace atlas { @@ -23,13 +25,13 @@ Mesh::Mesh(): Handle(new Implementation()) {} Mesh::Mesh(const Grid& grid, const eckit::Configuration& config): Handle([&]() { - if(config.has("mpi_comm")) { + if (config.has("mpi_comm")) { mpi::push(config.getString("mpi_comm")); } - util::Config cfg = grid.meshgenerator()|util::Config(config); - auto meshgenerator = MeshGenerator{grid.meshgenerator()|config}; - auto mesh = meshgenerator.generate(grid, grid::Partitioner(grid.partitioner()|config)); - if(config.has("mpi_comm")) { + auto cfg = grid.meshgenerator() | util::Config(config); + auto meshgenerator = MeshGenerator{grid.meshgenerator() | config}; + auto mesh = meshgenerator.generate(grid, grid::Partitioner(grid.partitioner() | config)); + if (config.has("mpi_comm")) { mpi::pop(); } mesh.get()->attach(); @@ -40,13 +42,13 @@ Mesh::Mesh(const Grid& grid, const eckit::Configuration& config): Mesh::Mesh(const Grid& grid, const grid::Partitioner& partitioner, const eckit::Configuration& config): Handle([&]() { - std::string mpi_comm = partitioner.mpi_comm(); - if(config.has("mpi_comm")) { + auto mpi_comm = partitioner.mpi_comm(); + if (config.has("mpi_comm")) { mpi_comm = config.getString("mpi_comm"); ATLAS_ASSERT(mpi_comm == partitioner.mpi_comm()); } mpi::Scope mpi_scope(mpi_comm); - auto meshgenerator = MeshGenerator{grid.meshgenerator()|config}; + auto meshgenerator = MeshGenerator{grid.meshgenerator() | config}; auto mesh = meshgenerator.generate(grid, partitioner); mesh.get()->attach(); return mesh.get(); @@ -57,7 +59,9 @@ Mesh::Mesh(const Grid& grid, const grid::Partitioner& partitioner, const eckit:: Mesh::Mesh(eckit::Stream& stream): Handle(new Implementation(stream)) {} -Mesh::operator bool() const { return get()->nodes().size() > 0; } +Mesh::operator bool() const { + return get()->nodes().size() > 0; +} //---------------------------------------------------------------------------------------------------------------------- diff --git a/src/atlas/projection/detail/MercatorProjection.cc b/src/atlas/projection/detail/MercatorProjection.cc index d3804523a..783ef74c8 100644 --- a/src/atlas/projection/detail/MercatorProjection.cc +++ b/src/atlas/projection/detail/MercatorProjection.cc @@ -10,6 +10,7 @@ #include #include +#include #include #include "eckit/config/Parametrisation.h" @@ -96,10 +97,18 @@ void MercatorProjectionT::lonlat2xy(double crd[]) const { rotation_.unrotate(crd); // then project - crd[XX] = k_radius_ * (D2R(normalise_mercator_(crd[LON]) - lon0_)); - crd[YY] = k_radius_ * 0.5 * std::log(t(crd[LAT])); - crd[XX] += false_easting_; - crd[YY] += false_northing_; + if (crd[LAT] >= 90. - 1e-3) { + crd[XX] = false_easting_; + crd[YY] = std::numeric_limits::infinity(); + } + else if (crd[LAT] <= -90. + 1e-3) { + crd[XX] = false_easting_; + crd[YY] = -std::numeric_limits::infinity(); + } + else { + crd[XX] = false_easting_ + k_radius_ * (D2R(normalise_mercator_(crd[LON]) - lon0_)); + crd[YY] = false_northing_ + k_radius_ * 0.5 * std::log(t(crd[LAT])); + } } template diff --git a/src/atlas/projection/detail/ProjectionImpl.cc b/src/atlas/projection/detail/ProjectionImpl.cc index 9b8264ac7..bf16eb28b 100644 --- a/src/atlas/projection/detail/ProjectionImpl.cc +++ b/src/atlas/projection/detail/ProjectionImpl.cc @@ -8,25 +8,23 @@ * nor does it submit to any jurisdiction. */ +#include "atlas/projection/detail/ProjectionImpl.h" + #include #include -#include +#include #include -#include +#include #include "eckit/types/FloatCompare.h" #include "eckit/utils/Hash.h" #include "eckit/utils/MD5.h" -#include "atlas/projection/detail/ProjectionImpl.h" - #include "atlas/projection/detail/ProjectionFactory.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" -namespace atlas { -namespace projection { -namespace detail { +namespace atlas::projection::detail { // -------------------------------------------------------------------------------------------------------------------- @@ -53,31 +51,19 @@ struct DerivateBuilder : public ProjectionImpl::DerivateFactory { struct DerivateForwards final : ProjectionImpl::Derivate { using Derivate::Derivate; - PointLonLat d(PointXY P) const override { - PointLonLat A(xy2lonlat(P)); - PointLonLat B(xy2lonlat(PointXY::add(P, H_))); - return PointLonLat::div(PointLonLat::sub(B, A), normH_); - } + PointLonLat d(PointXY P) const override { return (xy2lonlat(P + H_) - xy2lonlat(P)) * invnH_; } }; struct DerivateBackwards final : ProjectionImpl::Derivate { using Derivate::Derivate; - PointLonLat d(PointXY P) const override { - PointLonLat A(xy2lonlat(PointXY::sub(P, H_))); - PointLonLat B(xy2lonlat(P)); - return PointLonLat::div(PointLonLat::sub(B, A), normH_); - } + PointLonLat d(PointXY P) const override { return (xy2lonlat(P) - xy2lonlat(P - H_)) * invnH_; } }; struct DerivateCentral final : ProjectionImpl::Derivate { DerivateCentral(const ProjectionImpl& p, PointXY A, PointXY B, double h, double refLongitude): - Derivate(p, A, B, h, refLongitude), H2_{PointXY::mul(H_, 0.5)} {} + Derivate(p, A, B, h, refLongitude), H2_{H_ * 0.5} {} const PointXY H2_; - PointLonLat d(PointXY P) const override { - PointLonLat A(xy2lonlat(PointXY::sub(P, H2_))); - PointLonLat B(xy2lonlat(PointXY::add(P, H2_))); - return PointLonLat::div(PointLonLat::sub(B, A), normH_); - } + PointLonLat d(PointXY P) const override { return (xy2lonlat(P + H2_) - xy2lonlat(P - H2_)) * invnH_; } }; } // namespace @@ -85,7 +71,7 @@ struct DerivateCentral final : ProjectionImpl::Derivate { ProjectionImpl::Derivate::Derivate(const ProjectionImpl& p, PointXY A, PointXY B, double h, double refLongitude): projection_(p), H_{PointXY::mul(PointXY::normalize(PointXY::sub(B, A)), h)}, - normH_(PointXY::norm(H_)), + invnH_(1. / PointXY::norm(H_)), refLongitude_(refLongitude) {} ProjectionImpl::Derivate::~Derivate() = default; @@ -107,7 +93,7 @@ ProjectionImpl::Derivate* ProjectionImpl::DerivateFactory::build(const std::stri return new DerivateDegenerate(p, A, B, h, refLongitude); } - auto factory = get(type); + auto* factory = get(type); return factory->make(p, A, B, h); } @@ -184,7 +170,7 @@ ProjectionImpl::Normalise::Normalise(const eckit::Parametrisation& p) { provided = true; } if (provided) { - normalise_.reset(new util::NormaliseLongitude(values_[0], values_[1])); + normalise_ = std::make_unique(values_[0], values_[1]); } } @@ -192,7 +178,7 @@ ProjectionImpl::Normalise::Normalise(double west) { values_.resize(2); values_[0] = west; values_[1] = values_[0] + 360.; - normalise_.reset(new util::NormaliseLongitude(values_[0], values_[1])); + normalise_ = std::make_unique(values_[0], values_[1]); } @@ -246,8 +232,8 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) ATLAS_ASSERT(rect); // use central longitude as absolute reference (keep points within +-180 longitude range) - const auto centre = lonlat({(rect.xmin() + rect.xmax()) / 2., (rect.ymin() + rect.ymax()) / 2.}); - + const PointXY centre_xy{(rect.xmin() + rect.xmax()) / 2., (rect.ymin() + rect.ymax()) / 2.}; + const auto centre_lon = lonlat(centre_xy).lon(); const std::string derivative = "central"; constexpr double h_deg = 0.5e-6; // precision to microdegrees @@ -256,15 +242,18 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) const double h = units() == "degrees" ? h_deg : h_meters; + // 1. determine box from projected corners - const std::vector corners{ - {rect.xmin(), rect.ymax()}, {rect.xmax(), rect.ymax()}, {rect.xmax(), rect.ymin()}, {rect.xmin(), rect.ymin()}}; + const std::pair segments[] = {{{rect.xmin(), rect.ymax()}, {rect.xmax(), rect.ymax()}}, + {{rect.xmax(), rect.ymax()}, {rect.xmax(), rect.ymin()}}, + {{rect.xmax(), rect.ymin()}, {rect.xmin(), rect.ymin()}}, + {{rect.xmin(), rect.ymin()}, {rect.xmin(), rect.ymax()}}}; BoundLonLat bounds; - for (auto& p : corners) { + for (const auto [p, dummy] : segments) { auto q = lonlat(p); - longitude_in_range(centre.lon(), q.lon()); + longitude_in_range(centre_lon, q.lon()); bounds.extend(q, PointLonLat{h_deg, h_deg}); } @@ -272,12 +261,12 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) // 2. locate latitude extrema by checking if poles are included (in the un-projected frame) and if not, find extrema // not at the corners by refining iteratively - for (size_t i = 0; i < corners.size(); ++i) { - if (!bounds.includesNorthPole() || !bounds.includesSouthPole()) { - PointXY A = corners[i]; - PointXY B = corners[(i + 1) % corners.size()]; + bounds.includesNorthPole(bounds.includesNorthPole() || rect.contains(xy(PointLonLat{0, 90 - h_deg}))); + bounds.includesSouthPole(bounds.includesSouthPole() || rect.contains(xy(PointLonLat{0, -90 + h_deg}))); - std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre.lon())); + for (auto [A, B] : segments) { + if (!bounds.includesNorthPole() || !bounds.includesSouthPole()) { + std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre_lon)); double dAdy = derivate->d(A).lat(); double dBdy = derivate->d(B).lat(); @@ -307,12 +296,9 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) // 3. locate longitude extrema not at the corners by refining iteratively - for (size_t i = 0; i < corners.size(); ++i) { + for (auto [A, B] : segments) { if (!bounds.crossesDateLine()) { - PointXY A = corners[i]; - PointXY B = corners[(i + 1) % corners.size()]; - - std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre.lon())); + std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre_lon)); double dAdx = derivate->d(A).lon(); double dBdx = derivate->d(B).lon(); @@ -410,6 +396,4 @@ void atlas__Projection__lonlat2xy(const ProjectionImpl* This, const double lon, } // extern "C" -} // namespace detail -} // namespace projection -} // namespace atlas +} // namespace atlas::projection::detail diff --git a/src/atlas/projection/detail/ProjectionImpl.h b/src/atlas/projection/detail/ProjectionImpl.h index 6f777a5ee..559707643 100644 --- a/src/atlas/projection/detail/ProjectionImpl.h +++ b/src/atlas/projection/detail/ProjectionImpl.h @@ -10,13 +10,11 @@ #pragma once -#include #include #include #include #include "atlas/projection/Jacobian.h" -#include "atlas/runtime/Exception.h" #include "atlas/util/Factory.h" #include "atlas/util/NormaliseLongitude.h" #include "atlas/util/Object.h" @@ -36,9 +34,7 @@ class Config; } } // namespace atlas -namespace atlas { -namespace projection { -namespace detail { +namespace atlas::projection::detail { class ProjectionImpl : public util::Object { public: @@ -49,8 +45,7 @@ class ProjectionImpl : public util::Object { static const ProjectionImpl* create(const eckit::Parametrisation& p); static const ProjectionImpl* create(const std::string& type, const eckit::Parametrisation& p); - ProjectionImpl() = default; - virtual ~ProjectionImpl() = default; // destructor should be virtual + ProjectionImpl() = default; virtual std::string type() const = 0; @@ -123,7 +118,7 @@ class ProjectionImpl : public util::Object { protected: const ProjectionImpl& projection_; const PointXY H_; - const double normH_; + const double invnH_; const double refLongitude_; PointLonLat xy2lonlat(const PointXY& p) const; }; @@ -190,10 +185,8 @@ class NotRotated { static std::string classNamePrefix() { return ""; } // deliberately empty static std::string typePrefix() { return ""; } // deliberately empty - void rotate(double*) const { /* do nothing */ - } - void unrotate(double*) const { /* do nothing */ - } + void rotate(double*) const { /* do nothing */ } + void unrotate(double*) const { /* do nothing */ } bool rotated() const { return false; } @@ -211,6 +204,4 @@ void atlas__Projection__xy2lonlat(const ProjectionImpl* This, const double x, co void atlas__Projection__lonlat2xy(const ProjectionImpl* This, const double lon, const double lat, double& x, double& y); } -} // namespace detail -} // namespace projection -} // namespace atlas +} // namespace atlas::projection::detail From f1569ba16b05f9d63adedb527eeac60d56c2253a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 8 Apr 2024 17:49:03 +0200 Subject: [PATCH 22/23] Version 0.37.0 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 93d4c1ef0..0f1a7dfc7 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.36.0 +0.37.0 From d6f45e75c75d6868f9be9baebfd6a2203dcb7a07 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 9 Apr 2024 16:52:28 +0200 Subject: [PATCH 23/23] Update Changelog --- CHANGELOG.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 55f29a917..f7a0e8ece 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,19 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] +## [0.37.0] - 2024-04-09 +### Added +- Add SphericalVector interpolation method using parallel transport (#163) +- Added arrayForEachDim method + +### Fixed +- Bugfix for regional grids with ny > nx +- Fix build for configuration setting ATLAS_BITS_LOCAL=64 + +### Changed +- Use new LocalConfiguration baseclass functions in util::Config and util::Metadata instead of eckit::Value backdoor +- atlas_io is an adaptor library when eckit_codec is available (#181) + ## [0.36.0] - 2023-12-11 ### Added - Add TriangularMeshBuilder with Fortran API, so far for serial meshes only @@ -514,6 +527,7 @@ Fix StructuredInterpolation2D with retry for failed stencils ## 0.13.0 - 2018-02-16 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop +[0.37.0]: https://github.com/ecmwf/atlas/compare/0.36.0...0.37.0 [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