From e734e7a9a36417d0132e4a030069fd47bab6f0e0 Mon Sep 17 00:00:00 2001 From: Oliver Lomax Date: Fri, 15 Dec 2023 13:49:13 +0000 Subject: [PATCH] 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); }