From 2c9512fd9b122b0104042c2811a6042181b2b52f Mon Sep 17 00:00:00 2001 From: odlomax Date: Fri, 17 Nov 2023 16:29:04 +0000 Subject: [PATCH 01/26] Implemented forward interpolation. TODO: adjoint and more tests. --- src/atlas/CMakeLists.txt | 10 +- src/atlas/interpolation/method/Method.cc | 2 +- .../interpolation/method/MethodFactory.cc | 2 + .../paralleltransport/ParallelTransport.cc | 210 ++++++++++++++++++ .../paralleltransport/ParallelTransport.h | 62 ++++++ src/atlas/util/UnitSphere.h | 42 ++++ src/tests/interpolation/CMakeLists.txt | 8 +- .../test_interpolation_parallel_transport.cc | 171 ++++++++++++++ 8 files changed, 501 insertions(+), 6 deletions(-) create mode 100644 src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc create mode 100644 src/atlas/interpolation/method/paralleltransport/ParallelTransport.h create mode 100644 src/tests/interpolation/test_interpolation_parallel_transport.cc diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 6e0d3d01b..e68a4ad88 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -520,9 +520,11 @@ 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/method/paralleltransport/ParallelTransport.h +interpolation/method/paralleltransport/ParallelTransport.cc interpolation/Vector2D.cc interpolation/Vector2D.h interpolation/Vector3D.cc @@ -537,8 +539,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 @@ -862,7 +864,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/Method.cc b/src/atlas/interpolation/method/Method.cc index fb7f90221..2fc08f58c 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -374,7 +374,7 @@ void Method::do_execute(const FieldSet& fieldsSource, FieldSet& fieldsTarget, Me ATLAS_ASSERT(N == fieldsTarget.size()); for (idx_t i = 0; i < fieldsSource.size(); ++i) { - Method::do_execute(fieldsSource[i], fieldsTarget[i], metadata); + do_execute(fieldsSource[i], fieldsTarget[i], metadata); } } diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index 7c9f6952c..7a22e72c4 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -12,6 +12,7 @@ // for static linking #include "cubedsphere/CubedSphereBilinear.h" +#include "paralleltransport/ParallelTransport.h" #include "knn/GridBoxAverage.h" #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" @@ -47,6 +48,7 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); + MethodBuilder(); } } link; } diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc new file mode 100644 index 000000000..8a99603f0 --- /dev/null +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc @@ -0,0 +1,210 @@ +/* + * (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 +#include + +#include "atlas/array/ArrayView.h" +#include "atlas/array/helpers/ArrayForEach.h" +#include "atlas/array/Range.h" +#include "atlas/field/Field.h" +#include "atlas/interpolation/Cache.h" +#include "atlas/interpolation/Interpolation.h" +#include "atlas/interpolation/method/MethodFactory.h" +#include "atlas/interpolation/method/paralleltransport/ParallelTransport.h" +#include "atlas/linalg/sparse.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/UnitSphere.h" + +#include "eckit/linalg/Triplet.h" + + +namespace atlas { +namespace interpolation { +namespace method { + +namespace { +MethodBuilder __builder("parallel-transport"); + +template +void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { + + const auto nRows = matrix.rows(); + const auto nCols = matrix.cols(); + const auto rowIndices = matrix.outer(); + const auto colIndices = matrix.inner(); + auto valData = matrix.data(); + + atlas_omp_parallel_for (auto i = 0; i < nRows; ++i) { + for (auto dataIdx = rowIndices[i]; dataIdx < rowIndices[i+1]; ++dataIdx) { + const auto j = colIndices[dataIdx]; + auto&& value = valData[dataIdx]; + functor(value, i, j); + } + } +} + +} // namespace + +void ParallelTransport::do_setup(const Grid& source, const Grid& target, const Cache&) { + ATLAS_NOTIMPLEMENTED; +} + +void ParallelTransport::do_setup(const FunctionSpace& source, const FunctionSpace& target) { + ATLAS_TRACE("interpolation::method::ParallelTransport::do_setup"); + source_ = source; + target_ = target; + + if (target_.size() == 0) { + return; + } + + const auto baseInterpolator = Interpolation(interpolationScheme_, source_, target_); + setMatrix(MatrixCache(baseInterpolator)); + + // Get matrix dimensions. + const auto nRows = matrix().rows(); + const auto nCols = matrix().cols(); + + auto weights00 = std::vector(); + auto weights01 = std::vector(); + auto weights10 = std::vector(); + auto weights11 = std::vector(); + + const auto sourceLonLats = array::make_view(source_.lonlat()); + const auto targetLonLats = array::make_view(target_.lonlat()); + + + spaceMatrixForEach(matrix(), [&](auto&& weight, int i, int j){ + const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); + const auto targetLonLat = PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); + + const auto alpha = util::greatCircleCourse(sourceLonLat, targetLonLat); + + auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); + + + weights00.emplace_back(i, j, weight * std::cos(deltaAlpha)); + weights01.emplace_back(i, j, -weight * std::sin(deltaAlpha)); + weights10.emplace_back(i, j, weight * std::sin(deltaAlpha)); + weights11.emplace_back(i, j, weight * std::cos(deltaAlpha)); + }); + + + // Deal with slightly old fashioned Matrix interface + const auto buildMatrix = [&](auto& matrix, const auto& weights){ + auto tempMatrix = Matrix(nRows, nCols, weights); + matrix.swap(tempMatrix); + }; + + buildMatrix(matrix00_, weights00); + buildMatrix(matrix10_, weights10); + buildMatrix(matrix01_, weights01); + buildMatrix(matrix11_, weights11); + +} + +void ParallelTransport::print(std::ostream&) const { + ATLAS_NOTIMPLEMENTED; +} + +void ParallelTransport::do_execute(const Field &source, Field &target, Metadata &) const +{ + ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::do_execute()"); + + if (!(source.variables() == 2 || source.variables() == 3)) { + + auto metadata = Metadata(); + Method::do_execute(source, target, metadata); + + return; + + } + + if (target_.size() == 0) { + return; + } + + haloExchange(source); + + if (source.datatype().kind() == array::DataType::KIND_REAL64) { + interpolate_vector_field(source, target); + } + else if (source.datatype().kind() == array::DataType::KIND_REAL32) { + interpolate_vector_field(source, target); + } + else { + ATLAS_NOTIMPLEMENTED; + } + + target.set_dirty(); +} + +template +void ParallelTransport::interpolate_vector_field(const Field &source, Field &target) const +{ + if (source.rank() == 2) { + interpolate_vector_field(source, target); + } + else if (source.rank() == 3 ) { + interpolate_vector_field(source, target); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + +template +void ParallelTransport::interpolate_vector_field(const Field &source, Field &target) const +{ + using namespace linalg; + + const auto sourceView = array::make_view(source); + auto targetView = array::make_view(target); + + array::make_view(target).assign(0); + + + + const auto matrixMultiply = [&](const auto& matrix, int sourceVariableIdx, int targetVariableIdx) { + + spaceMatrixForEach(matrix, [&](const auto& weight, int i, int j){ + + const auto adder = [&](auto&& targetElem, auto&& sourceElem){ + targetElem += weight * sourceElem; + }; + + if constexpr (Rank == 2) { + adder(targetView(i, targetVariableIdx), sourceView(j, sourceVariableIdx)); + } + else { + const auto sourceSlice = sourceView.slice(j, array::Range::all(), sourceVariableIdx); + auto targetSlice = targetView.slice(i, array::Range::all(), targetVariableIdx); + array::helpers::ArrayForEach<0>::apply(std::tie(targetSlice, sourceSlice), adder); + } + + }); + }; + + matrixMultiply(matrix00_, 0, 0); + matrixMultiply(matrix01_, 1, 0); + matrixMultiply(matrix10_, 0, 1); + matrixMultiply(matrix11_, 1, 1); + + if (source.variables() == 3) { + matrixMultiply(matrix(), 2, 2); + } + +} + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h new file mode 100644 index 000000000..3a71db20d --- /dev/null +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h @@ -0,0 +1,62 @@ +/* + * (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. + */ + +#pragma once + +#include "atlas/functionspace/FunctionSpace.h" +#include "atlas/interpolation/method/Method.h" +#include "atlas/linalg/sparse.h" +#include "eckit/config/Configuration.h" + +namespace atlas { +namespace interpolation { +namespace method { + +class ParallelTransport : public Method { + public: + ParallelTransport(const Config& config): Method(config) { + const auto& conf = dynamic_cast(config); + interpolationScheme_ = conf.getSubConfiguration("scheme"); + + } + virtual ~ParallelTransport() 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 Field& source, Field& target, 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_; + + // Matrices for vector field interpolation + Matrix matrix00_; + Matrix matrix01_; + Matrix matrix10_; + Matrix matrix11_; + +}; + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/util/UnitSphere.h b/src/atlas/util/UnitSphere.h index c1458468f..1f38e455d 100644 --- a/src/atlas/util/UnitSphere.h +++ b/src/atlas/util/UnitSphere.h @@ -10,6 +10,12 @@ #pragma once +#include +#include + +#include "atlas/util/Constants.h" +#include "atlas/util/Point.h" + #include "eckit/geometry/UnitSphere.h" //------------------------------------------------------------------------------------------------------ @@ -21,6 +27,42 @@ namespace util { using eckit::geometry::UnitSphere; +/// @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 +/// +inline std::pair greatCircleCourse(const Point2& lonLat1, + const Point2& lonLat2) { + + const auto lambda1 = lonLat1[0] * Constants::degreesToRadians(); + const auto lambda2 = lonLat2[0] * Constants::degreesToRadians(); + const auto phi1 = lonLat1[1] * Constants::degreesToRadians(); + const auto phi2 = lonLat2[1] * 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 * Constants::radiansToDegrees(), + alpha2 * Constants::radiansToDegrees()); +}; + //------------------------------------------------------------------------------------------------------ } // namespace util diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 596729209..4e7edf257 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -81,6 +81,12 @@ ecbuild_add_test( TARGET atlas_test_interpolation_non_linear ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +ecbuild_add_test( TARGET atlas_test_interpolation_parallel_transport + SOURCES test_interpolation_parallel_transport.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + ecbuild_add_test( TARGET atlas_test_interpolation_bilinear COMMAND atlas_test_interpolation_structured2D ARGS --scheme linear ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} @@ -122,4 +128,4 @@ ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -endif() \ No newline at end of file +endif() diff --git a/src/tests/interpolation/test_interpolation_parallel_transport.cc b/src/tests/interpolation/test_interpolation_parallel_transport.cc new file mode 100644 index 000000000..65029f5bf --- /dev/null +++ b/src/tests/interpolation/test_interpolation_parallel_transport.cc @@ -0,0 +1,171 @@ +/* + * (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/function/VortexRollup.h" +#include "atlas/util/UnitSphere.h" + +#include "tests/AtlasTestEnvironment.h" + +namespace atlas { +namespace test { + +using namespace atlas::util; +using namespace atlas::array::helpers; + +// 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) { + + // set hLon and hLat step size. + const double hLon = 0.0001; + const double hLat = 0.0001; + + // Get finite differences. + + // Set u. + const double u = (util::function::vortex_rollup(lon, lat + 0.5 * hLat, 0.1) - + util::function::vortex_rollup(lon, lat - 0.5 * hLat, 0.1)) / + hLat; + + const double v = -(util::function::vortex_rollup(lon + 0.5 * hLon, lat, 0.1) - + util::function::vortex_rollup(lon - 0.5 * hLon, lat, 0.1)) / + (hLon * std::cos(lat * util::Constants::degreesToRadians())); + + return std::make_pair(u, v); +} + +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); +} + +void testInterpolation(const Config& config) { + + const auto sourceGrid = Grid(config.getString("source_grid")); + const auto sourceMesh = + MeshGenerator(config.getString("source_mesh")).generate(sourceGrid); + const auto sourceFunctionSpace = functionspace::NodeColumns(sourceMesh); + + const auto targetGrid = Grid(config.getString("target_grid")); + const auto targetMesh = + MeshGenerator(config.getString("target_mesh")).generate(targetGrid); + const auto targetFunctionSpace = functionspace::NodeColumns(targetMesh); + + auto sourceFieldSet = FieldSet{}; + auto targetFieldSet = FieldSet{}; + + const auto sourceLonLat = + array::make_view(sourceFunctionSpace.lonlat()); + const auto targetLonLat = + array::make_view(targetFunctionSpace.lonlat()); + + auto sourceField = array::make_view(sourceFieldSet.add( + sourceFunctionSpace.createField(option::name("test field") | + option::levels(1) | + option::variables(2)))); + + auto targetField = array::make_view(targetFieldSet.add( + targetFunctionSpace.createField(option::name("test field") | + option::levels(1) | + option::variables(2)))); + + ArrayForEach<0>::apply(std::tie(sourceLonLat, sourceField), + [](auto&& lonLat, auto&& sourceColumn) { + ArrayForEach<0>::apply(std::tie(sourceColumn), [&](auto&& sourceElem) { + std::tie(sourceElem(0), sourceElem(1)) = + vortexField(lonLat(0), lonLat(1)); + }); + }); + + const auto interp = Interpolation(config.getSubConfiguration("scheme"), + sourceFunctionSpace, targetFunctionSpace); + + interp.execute(sourceFieldSet, targetFieldSet); + targetFieldSet.haloExchange(); + + auto errorField = array::make_view( + targetFieldSet.add(targetFunctionSpace.createField( + option::name("error field") | option::levels(1)))); + + ArrayForEach<0>::apply( + std::tie(targetLonLat, targetField, errorField), + [](auto&& lonLat, auto&& targetColumn, auto&& errorColumn) { + ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), + [&](auto&& targetElem, auto&& errorElem) { + + auto deltaVal = vortexField(lonLat(0), lonLat(1)); + deltaVal.first -= targetElem(0); + deltaVal.second -= targetElem(1); + + errorElem = std::sqrt(deltaVal.first * deltaVal.first + + deltaVal.second * deltaVal.second); + }); + }); + + gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); + gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); +} + +CASE("great-circle course") { + + 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 ] = greatCircleCourse(point1, point2); + const auto[ course3, course4 ] = greatCircleCourse(point3, point4); + + EXPECT_APPROX_EQ(course1, targetCourse1, 0.01); + EXPECT_APPROX_EQ(course2, targetCourse2, 0.01); + EXPECT_APPROX_EQ(course3, targetCourse3, 1.e-7); + EXPECT_APPROX_EQ(std::abs(course4), targetCourse4, 1.e-7); +} + +CASE("cubed sphere vector interpolation") { + + const auto baseInterpScheme = + util::Config("type", "cubedsphere-bilinear").set("adjoint", true); + const auto interpScheme = util::Config("type", "parallel-transport") + .set("scheme", baseInterpScheme); + const auto cubedSphereConf = Config("source_grid", "CS-LFR-48") + .set("source_mesh", "cubedsphere_dual") + .set("target_grid", "O48") + .set("target_mesh", "structured") + .set("file_id", "parallel_transport_cs") + .set("scheme", baseInterpScheme); + + testInterpolation((cubedSphereConf)); +} +} +} + +int main(int argc, char** argv) { return atlas::test::run(argc, argv); } From 7c0a5d7024c04e27fa1376ab6b484f5021671270 Mon Sep 17 00:00:00 2001 From: odlomax Date: Mon, 20 Nov 2023 18:36:48 +0000 Subject: [PATCH 02/26] Fixed race condition errors. --- src/atlas/interpolation/method/Method.cc | 2 +- .../paralleltransport/ParallelTransport.cc | 97 ++++++++++++------- .../paralleltransport/ParallelTransport.h | 3 +- src/tests/interpolation/CMakeLists.txt | 1 + .../test_interpolation_parallel_transport.cc | 17 ++++ 5 files changed, 82 insertions(+), 38 deletions(-) diff --git a/src/atlas/interpolation/method/Method.cc b/src/atlas/interpolation/method/Method.cc index 2fc08f58c..25701dbe6 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -374,7 +374,7 @@ void Method::do_execute(const FieldSet& fieldsSource, FieldSet& fieldsTarget, Me ATLAS_ASSERT(N == fieldsTarget.size()); for (idx_t i = 0; i < fieldsSource.size(); ++i) { - do_execute(fieldsSource[i], fieldsTarget[i], metadata); + Method::do_execute(fieldsSource[i], fieldsTarget[i], metadata); } } diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc index 8a99603f0..30c3820f3 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc @@ -12,6 +12,7 @@ #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" @@ -43,11 +44,20 @@ void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { const auto colIndices = matrix.inner(); auto valData = matrix.data(); - atlas_omp_parallel_for (auto i = 0; i < nRows; ++i) { + atlas_omp_parallel_for (auto i = size_t{}; i < nRows; ++i) { for (auto dataIdx = rowIndices[i]; dataIdx < rowIndices[i+1]; ++dataIdx) { - const auto j = colIndices[dataIdx]; - auto&& value = valData[dataIdx]; - functor(value, i, j); + const auto j = size_t(colIndices[dataIdx]); + auto& value = valData[dataIdx]; + + if constexpr (std::is_invocable_v) { + functor(value, i, j); + } + else if constexpr (std::is_invocable_v) { + functor(value, i, j, dataIdx); + } + else { + ATLAS_NOTIMPLEMENTED; + } } } } @@ -73,17 +83,18 @@ void ParallelTransport::do_setup(const FunctionSpace& source, const FunctionSpac // Get matrix dimensions. const auto nRows = matrix().rows(); const auto nCols = matrix().cols(); + const auto nNonZeros = matrix().nonZeros(); - auto weights00 = std::vector(); - auto weights01 = std::vector(); - auto weights10 = std::vector(); - auto weights11 = std::vector(); + auto weights00 = std::vector(nNonZeros); + auto weights01 = std::vector(nNonZeros); + auto weights10 = std::vector(nNonZeros); + auto weights11 = std::vector(nNonZeros); const auto sourceLonLats = array::make_view(source_.lonlat()); const auto targetLonLats = array::make_view(target_.lonlat()); - spaceMatrixForEach(matrix(), [&](auto&& weight, int i, int j){ + spaceMatrixForEach(matrix(), [&](auto&& weight, auto i, auto j, auto dataIdx){ const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); const auto targetLonLat = PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); @@ -91,11 +102,10 @@ void ParallelTransport::do_setup(const FunctionSpace& source, const FunctionSpac auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - - weights00.emplace_back(i, j, weight * std::cos(deltaAlpha)); - weights01.emplace_back(i, j, -weight * std::sin(deltaAlpha)); - weights10.emplace_back(i, j, weight * std::sin(deltaAlpha)); - weights11.emplace_back(i, j, weight * std::cos(deltaAlpha)); + weights00[dataIdx] = {i, j, weight * std::cos(deltaAlpha)}; + weights01[dataIdx] = {i, j, -weight * std::sin(deltaAlpha)}; + weights10[dataIdx] = {i, j, weight * std::sin(deltaAlpha)}; + weights11[dataIdx] = {i, j, weight * std::cos(deltaAlpha)}; }); @@ -116,14 +126,26 @@ void ParallelTransport::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } -void ParallelTransport::do_execute(const Field &source, Field &target, Metadata &) const +void ParallelTransport::do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, Metadata& metadata) const +{ + ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::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 ParallelTransport::do_execute(const Field& sourceField, Field& targetField, Metadata &) const { ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::do_execute()"); - if (!(source.variables() == 2 || source.variables() == 3)) { + if (!(sourceField.variables() == 2 || sourceField.variables() == 3)) { auto metadata = Metadata(); - Method::do_execute(source, target, metadata); + Method::do_execute(sourceField, targetField, metadata); return; @@ -133,29 +155,29 @@ void ParallelTransport::do_execute(const Field &source, Field &target, Metadata return; } - haloExchange(source); + haloExchange(sourceField); - if (source.datatype().kind() == array::DataType::KIND_REAL64) { - interpolate_vector_field(source, target); + if (sourceField.datatype().kind() == array::DataType::KIND_REAL64) { + interpolate_vector_field(sourceField, targetField); } - else if (source.datatype().kind() == array::DataType::KIND_REAL32) { - interpolate_vector_field(source, target); + else if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { + interpolate_vector_field(sourceField, targetField); } else { ATLAS_NOTIMPLEMENTED; } - target.set_dirty(); + targetField.set_dirty(); } template -void ParallelTransport::interpolate_vector_field(const Field &source, Field &target) const +void ParallelTransport::interpolate_vector_field(const Field& sourceField, Field& targetField) const { - if (source.rank() == 2) { - interpolate_vector_field(source, target); + if (sourceField.rank() == 2) { + interpolate_vector_field(sourceField, targetField); } - else if (source.rank() == 3 ) { - interpolate_vector_field(source, target); + else if (sourceField.rank() == 3 ) { + interpolate_vector_field(sourceField, targetField); } else { ATLAS_NOTIMPLEMENTED; @@ -163,20 +185,20 @@ void ParallelTransport::interpolate_vector_field(const Field &source, Field &tar } template -void ParallelTransport::interpolate_vector_field(const Field &source, Field &target) const +void ParallelTransport::interpolate_vector_field(const Field& sourceField, Field& targetField) const { using namespace linalg; - const auto sourceView = array::make_view(source); - auto targetView = array::make_view(target); + const auto sourceView = array::make_view(sourceField); + auto targetView = array::make_view(targetField); - array::make_view(target).assign(0); + array::make_view(targetField).assign(0); - const auto matrixMultiply = [&](const auto& matrix, int sourceVariableIdx, int targetVariableIdx) { + const auto matrixMultiply = [&](const auto& matrix, auto sourceVariableIdx, auto targetVariableIdx) { - spaceMatrixForEach(matrix, [&](const auto& weight, int i, int j){ + spaceMatrixForEach(matrix, [&](const auto& weight, auto i, auto j){ const auto adder = [&](auto&& targetElem, auto&& sourceElem){ targetElem += weight * sourceElem; @@ -185,11 +207,14 @@ void ParallelTransport::interpolate_vector_field(const Field &source, Field &tar if constexpr (Rank == 2) { adder(targetView(i, targetVariableIdx), sourceView(j, sourceVariableIdx)); } - else { + else if constexpr (Rank == 3) { const auto sourceSlice = sourceView.slice(j, array::Range::all(), sourceVariableIdx); auto targetSlice = targetView.slice(i, array::Range::all(), targetVariableIdx); array::helpers::ArrayForEach<0>::apply(std::tie(targetSlice, sourceSlice), adder); } + else { + ATLAS_NOTIMPLEMENTED; + } }); }; @@ -199,7 +224,7 @@ void ParallelTransport::interpolate_vector_field(const Field &source, Field &tar matrixMultiply(matrix10_, 0, 1); matrixMultiply(matrix11_, 1, 1); - if (source.variables() == 3) { + if (sourceField.variables() == 3) { matrixMultiply(matrix(), 2, 2); } diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h index 3a71db20d..7384077c8 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h @@ -29,7 +29,8 @@ class ParallelTransport : public Method { const FunctionSpace& source() const override { return source_; } const FunctionSpace& target() const override { return target_; } - void do_execute(const Field& source, Field& target, Metadata&) const override; + 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: diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 4e7edf257..c680db212 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -82,6 +82,7 @@ ecbuild_add_test( TARGET atlas_test_interpolation_non_linear ) ecbuild_add_test( TARGET atlas_test_interpolation_parallel_transport + OMP 4 SOURCES test_interpolation_parallel_transport.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} diff --git a/src/tests/interpolation/test_interpolation_parallel_transport.cc b/src/tests/interpolation/test_interpolation_parallel_transport.cc index 65029f5bf..98b87dcea 100644 --- a/src/tests/interpolation/test_interpolation_parallel_transport.cc +++ b/src/tests/interpolation/test_interpolation_parallel_transport.cc @@ -161,10 +161,27 @@ CASE("cubed sphere vector interpolation") { .set("target_grid", "O48") .set("target_mesh", "structured") .set("file_id", "parallel_transport_cs") + .set("scheme", interpScheme); + + testInterpolation((cubedSphereConf)); +} + +CASE("finite element vector interpolation") { + + const auto baseInterpScheme = + util::Config("type", "finite-element").set("adjoint", true); + const auto interpScheme = util::Config("type", "parallel-transport") + .set("scheme", baseInterpScheme); + const auto cubedSphereConf = Config("source_grid", "O48") + .set("source_mesh", "structured") + .set("target_grid", "CS-LFR-48") + .set("target_mesh", "cubedsphere_dual") + .set("file_id", "parallel_transport_fe") .set("scheme", baseInterpScheme); testInterpolation((cubedSphereConf)); } + } } From d01f5dd07cfb8831fbe179d217f6ebf43ed2fa7b Mon Sep 17 00:00:00 2001 From: odlomax Date: Tue, 21 Nov 2023 12:15:27 +0000 Subject: [PATCH 03/26] Replaced rotation matrices with complex interpolation weights. --- .../paralleltransport/ParallelTransport.cc | 197 +++++++++--------- .../paralleltransport/ParallelTransport.h | 9 +- .../test_interpolation_parallel_transport.cc | 2 +- 3 files changed, 109 insertions(+), 99 deletions(-) diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc index 30c3820f3..12f68b300 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc @@ -27,7 +27,6 @@ #include "eckit/linalg/Triplet.h" - namespace atlas { namespace interpolation { namespace method { @@ -35,7 +34,7 @@ namespace method { namespace { MethodBuilder __builder("parallel-transport"); -template +template void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { const auto nRows = matrix.rows(); @@ -44,17 +43,21 @@ void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { const auto colIndices = matrix.inner(); auto valData = matrix.data(); - atlas_omp_parallel_for (auto i = size_t{}; i < nRows; ++i) { - for (auto dataIdx = rowIndices[i]; dataIdx < rowIndices[i+1]; ++dataIdx) { + atlas_omp_parallel_for(auto i = size_t{}; i < nRows; ++i) { + for (auto dataIdx = rowIndices[i]; dataIdx < rowIndices[i + 1]; ++dataIdx) { const auto j = size_t(colIndices[dataIdx]); - auto& value = valData[dataIdx]; - - if constexpr (std::is_invocable_v) { - functor(value, i, j); - } - else if constexpr (std::is_invocable_v) { - functor(value, i, j, dataIdx); - } + auto&& value = valData[dataIdx]; + + if + constexpr( + std::is_invocable_v) { + functor(value, i, j); + } + else if + constexpr(std::is_invocable_v) { + functor(value, i, j, dataIdx); + } else { ATLAS_NOTIMPLEMENTED; } @@ -62,13 +65,47 @@ void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { } } +template +void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, + TargetView&& targetView, const Functor& mappingFunctor) { + + spaceMatrixForEach(matrix, [&](const auto& weight, auto i, auto j) { + + 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()); + mappingFunctor(weight, sourceSlice, targetSlice); + } + else if + constexpr(rank == 3) { + const auto iterationFuctor = [&](auto&& sourceVars, auto&& targetVars) { + mappingFunctor(weight, sourceVars, targetVars); + }; + 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), iterationFuctor); + } + else { + ATLAS_NOTIMPLEMENTED; + } + }); +} + } // namespace -void ParallelTransport::do_setup(const Grid& source, const Grid& target, const Cache&) { +void ParallelTransport::do_setup(const Grid& source, const Grid& target, + const Cache&) { ATLAS_NOTIMPLEMENTED; } -void ParallelTransport::do_setup(const FunctionSpace& source, const FunctionSpace& target) { +void ParallelTransport::do_setup(const FunctionSpace& source, + const FunctionSpace& target) { ATLAS_TRACE("interpolation::method::ParallelTransport::do_setup"); source_ = source; target_ = target; @@ -77,7 +114,8 @@ void ParallelTransport::do_setup(const FunctionSpace& source, const FunctionSpac return; } - const auto baseInterpolator = Interpolation(interpolationScheme_, source_, target_); + const auto baseInterpolator = + Interpolation(interpolationScheme_, source_, target_); setMatrix(MatrixCache(baseInterpolator)); // Get matrix dimensions. @@ -85,49 +123,44 @@ void ParallelTransport::do_setup(const FunctionSpace& source, const FunctionSpac const auto nCols = matrix().cols(); const auto nNonZeros = matrix().nonZeros(); - auto weights00 = std::vector(nNonZeros); - auto weights01 = std::vector(nNonZeros); - auto weights10 = std::vector(nNonZeros); - auto weights11 = std::vector(nNonZeros); + auto weightsReal = std::vector(nNonZeros); + auto weightsImag = std::vector(nNonZeros); const auto sourceLonLats = array::make_view(source_.lonlat()); const auto targetLonLats = array::make_view(target_.lonlat()); - - spaceMatrixForEach(matrix(), [&](auto&& weight, auto i, auto j, auto dataIdx){ - const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); - const auto targetLonLat = PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); + // Make complex weights (would be nice if we could have a complex matrix). + spaceMatrixForEach(matrix(), + [&](auto&& weight, auto i, auto j, auto dataIdx) { + const auto sourceLonLat = + PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); + const auto targetLonLat = + PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); const auto alpha = util::greatCircleCourse(sourceLonLat, targetLonLat); - auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); + auto deltaAlpha = + (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - weights00[dataIdx] = {i, j, weight * std::cos(deltaAlpha)}; - weights01[dataIdx] = {i, j, -weight * std::sin(deltaAlpha)}; - weights10[dataIdx] = {i, j, weight * std::sin(deltaAlpha)}; - weights11[dataIdx] = {i, j, weight * std::cos(deltaAlpha)}; + weightsReal[dataIdx] = {i, j, weight * std::cos(deltaAlpha)}; + weightsImag[dataIdx] = {i, j, weight * std::sin(deltaAlpha)}; }); - // Deal with slightly old fashioned Matrix interface - const auto buildMatrix = [&](auto& matrix, const auto& weights){ + const auto buildMatrix = [&](auto& matrix, const auto& weights) { auto tempMatrix = Matrix(nRows, nCols, weights); matrix.swap(tempMatrix); }; - buildMatrix(matrix00_, weights00); - buildMatrix(matrix10_, weights10); - buildMatrix(matrix01_, weights01); - buildMatrix(matrix11_, weights11); - + buildMatrix(matrixReal_, weightsReal); + buildMatrix(matrixImag_, weightsImag); } -void ParallelTransport::print(std::ostream&) const { - ATLAS_NOTIMPLEMENTED; -} +void ParallelTransport::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } -void ParallelTransport::do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, Metadata& metadata) const -{ +void ParallelTransport::do_execute(const FieldSet& sourceFieldSet, + FieldSet& targetFieldSet, + Metadata& metadata) const { ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::do_execute()"); const auto nFields = sourceFieldSet.size(); @@ -138,8 +171,8 @@ void ParallelTransport::do_execute(const FieldSet& sourceFieldSet, FieldSet& tar } } -void ParallelTransport::do_execute(const Field& sourceField, Field& targetField, Metadata &) const -{ +void ParallelTransport::do_execute(const Field& sourceField, Field& targetField, + Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::do_execute()"); if (!(sourceField.variables() == 2 || sourceField.variables() == 3)) { @@ -148,7 +181,6 @@ void ParallelTransport::do_execute(const Field& sourceField, Field& targetField, Method::do_execute(sourceField, targetField, metadata); return; - } if (target_.size() == 0) { @@ -159,75 +191,54 @@ void ParallelTransport::do_execute(const Field& sourceField, Field& targetField, if (sourceField.datatype().kind() == array::DataType::KIND_REAL64) { interpolate_vector_field(sourceField, targetField); - } - else if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { + } else if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { interpolate_vector_field(sourceField, targetField); - } - else { + } else { ATLAS_NOTIMPLEMENTED; } targetField.set_dirty(); } -template -void ParallelTransport::interpolate_vector_field(const Field& sourceField, Field& targetField) const -{ +template +void ParallelTransport::interpolate_vector_field(const Field& sourceField, + Field& targetField) const { if (sourceField.rank() == 2) { - interpolate_vector_field(sourceField, targetField); - } - else if (sourceField.rank() == 3 ) { + interpolate_vector_field(sourceField, targetField); + } else if (sourceField.rank() == 3) { interpolate_vector_field(sourceField, targetField); - } - else { + } else { ATLAS_NOTIMPLEMENTED; } } -template -void ParallelTransport::interpolate_vector_field(const Field& sourceField, Field& targetField) const -{ - using namespace linalg; +template +void ParallelTransport::interpolate_vector_field(const Field& sourceField, + Field& targetField) const { const auto sourceView = array::make_view(sourceField); auto targetView = array::make_view(targetField); - - array::make_view(targetField).assign(0); - - - - const auto matrixMultiply = [&](const auto& matrix, auto sourceVariableIdx, auto targetVariableIdx) { - - spaceMatrixForEach(matrix, [&](const auto& weight, auto i, auto j){ - - const auto adder = [&](auto&& targetElem, auto&& sourceElem){ - targetElem += weight * sourceElem; - }; - - if constexpr (Rank == 2) { - adder(targetView(i, targetVariableIdx), sourceView(j, sourceVariableIdx)); - } - else if constexpr (Rank == 3) { - const auto sourceSlice = sourceView.slice(j, array::Range::all(), sourceVariableIdx); - auto targetSlice = targetView.slice(i, array::Range::all(), targetVariableIdx); - array::helpers::ArrayForEach<0>::apply(std::tie(targetSlice, sourceSlice), adder); - } - else { - ATLAS_NOTIMPLEMENTED; - } - - }); - }; - - matrixMultiply(matrix00_, 0, 0); - matrixMultiply(matrix01_, 1, 0); - matrixMultiply(matrix10_, 0, 1); - matrixMultiply(matrix11_, 1, 1); + targetView.assign(0.); + + // Matrix multiplication split in two to simulate complex variable + // multiplication. + matrixMultiply(matrixReal_, sourceView, targetView, + [](const auto& weight, auto&& sourceVars, auto&& targetVars) { + targetVars(0) += weight * sourceVars(0); + targetVars(1) += weight * sourceVars(1); + }); + matrixMultiply(matrixImag_, sourceView, targetView, + [](const auto& weight, auto&& sourceVars, auto&& targetVars) { + targetVars(0) -= weight * sourceVars(1); + targetVars(1) += weight * sourceVars(0); + }); if (sourceField.variables() == 3) { - matrixMultiply(matrix(), 2, 2); + matrixMultiply( + matrix(), sourceView, targetView, + [](const auto& weight, auto&& sourceVars, + auto&& targetVars) { targetVars(2) = weight * sourceVars(2); }); } - } } // namespace method diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h index 7384077c8..fbc427104 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h @@ -50,11 +50,10 @@ class ParallelTransport : public Method { FunctionSpace source_; FunctionSpace target_; - // Matrices for vector field interpolation - Matrix matrix00_; - Matrix matrix01_; - Matrix matrix10_; - Matrix matrix11_; + // Complex interpolation weights. We treat a (u, v) pair as a complex variable + // with u on the real number line and v on the imaginary number line. + Matrix matrixReal_; + Matrix matrixImag_; }; diff --git a/src/tests/interpolation/test_interpolation_parallel_transport.cc b/src/tests/interpolation/test_interpolation_parallel_transport.cc index 98b87dcea..23967c6f8 100644 --- a/src/tests/interpolation/test_interpolation_parallel_transport.cc +++ b/src/tests/interpolation/test_interpolation_parallel_transport.cc @@ -177,7 +177,7 @@ CASE("finite element vector interpolation") { .set("target_grid", "CS-LFR-48") .set("target_mesh", "cubedsphere_dual") .set("file_id", "parallel_transport_fe") - .set("scheme", baseInterpScheme); + .set("scheme", interpScheme); testInterpolation((cubedSphereConf)); } From d5c9ef5db656add22fc522a313ebe6ca00ea2383 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 23 Nov 2023 11:30:29 +0000 Subject: [PATCH 04/26] Fixed typos before discussion. --- .../paralleltransport/ParallelTransport.cc | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc index 12f68b300..59e23c18e 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc +++ b/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc @@ -35,7 +35,7 @@ namespace { MethodBuilder __builder("parallel-transport"); template -void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { +void sparseMatrixForEach(MatrixT&& matrix, const Functor& functor) { const auto nRows = matrix.rows(); const auto nCols = matrix.cols(); @@ -48,14 +48,12 @@ void spaceMatrixForEach(MatrixT&& matrix, const Functor& functor) { const auto j = size_t(colIndices[dataIdx]); auto&& value = valData[dataIdx]; - if - constexpr( + if constexpr( std::is_invocable_v) { functor(value, i, j); } - else if - constexpr(std::is_invocable_v) { + else if constexpr(std::is_invocable_v) { functor(value, i, j, dataIdx); } else { @@ -70,17 +68,15 @@ template ::rank(); - if - constexpr(rank == 2) { + if constexpr(rank == 2) { const auto sourceSlice = sourceView.slice(j, array::Range::all()); auto targetSlice = targetView.slice(i, array::Range::all()); mappingFunctor(weight, sourceSlice, targetSlice); } - else if - constexpr(rank == 3) { + else if constexpr(rank == 3) { const auto iterationFuctor = [&](auto&& sourceVars, auto&& targetVars) { mappingFunctor(weight, sourceVars, targetVars); }; @@ -130,7 +126,7 @@ void ParallelTransport::do_setup(const FunctionSpace& source, const auto targetLonLats = array::make_view(target_.lonlat()); // Make complex weights (would be nice if we could have a complex matrix). - spaceMatrixForEach(matrix(), + sparseMatrixForEach(matrix(), [&](auto&& weight, auto i, auto j, auto dataIdx) { const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); From 088403efd3bcd439bc0558cc311a87857e5ed81b Mon Sep 17 00:00:00 2001 From: odlomax Date: Tue, 28 Nov 2023 17:33:31 +0000 Subject: [PATCH 05/26] Renamed class to SphericalVector. Added Eigen3 matrices. --- src/atlas/CMakeLists.txt | 11 +- .../interpolation/method/MethodFactory.cc | 4 +- .../SphericalVector.cc} | 130 ++++++++---------- .../SphericalVector.h} | 22 ++- src/tests/interpolation/CMakeLists.txt | 5 +- ...=> test_interpolation_spherical_vector.cc} | 8 +- 6 files changed, 88 insertions(+), 92 deletions(-) rename src/atlas/interpolation/method/{paralleltransport/ParallelTransport.cc => sphericalvector/SphericalVector.cc} (56%) rename src/atlas/interpolation/method/{paralleltransport/ParallelTransport.h => sphericalvector/SphericalVector.h} (73%) rename src/tests/interpolation/{test_interpolation_parallel_transport.cc => test_interpolation_spherical_vector.cc} (95%) diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 366fdc1e5..4a18dec2a 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -525,8 +525,6 @@ functionspace/detail/CubedSphereStructure.cc # for cubedsphere matching mesh partitioner interpolation/method/cubedsphere/CellFinder.cc interpolation/method/cubedsphere/CellFinder.h -interpolation/method/paralleltransport/ParallelTransport.h -interpolation/method/paralleltransport/ParallelTransport.cc interpolation/Vector2D.cc interpolation/Vector2D.h interpolation/Vector3D.cc @@ -671,6 +669,15 @@ interpolation/nonlinear/NonLinear.cc interpolation/nonlinear/NonLinear.h ) +if (eckit_EIGEN_FOUND) + +list (APPEND atlas_interpolation_srcs +interpolation/method/sphericalvector/SphericalVector.h +interpolation/method/sphericalvector/SphericalVector.cc +) + +endif() + list( APPEND atlas_linalg_srcs linalg/Indexing.h linalg/Introspection.h diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index 7a22e72c4..9f2bff494 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -12,7 +12,7 @@ // for static linking #include "cubedsphere/CubedSphereBilinear.h" -#include "paralleltransport/ParallelTransport.h" +#include "sphericalvector/SphericalVector.h" #include "knn/GridBoxAverage.h" #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" @@ -48,7 +48,7 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); - MethodBuilder(); + MethodBuilder(); } } link; } diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc similarity index 56% rename from src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc rename to src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 59e23c18e..4376a98ae 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -16,7 +16,7 @@ #include "atlas/interpolation/Cache.h" #include "atlas/interpolation/Interpolation.h" #include "atlas/interpolation/method/MethodFactory.h" -#include "atlas/interpolation/method/paralleltransport/ParallelTransport.h" +#include "atlas/interpolation/method/sphericalvector/SphericalVector.h" #include "atlas/linalg/sparse.h" #include "atlas/option/Options.h" #include "atlas/parallel/omp/omp.h" @@ -32,33 +32,14 @@ namespace interpolation { namespace method { namespace { -MethodBuilder __builder("parallel-transport"); +MethodBuilder __builder("spherical-vector"); template -void sparseMatrixForEach(MatrixT&& matrix, const Functor& functor) { - - const auto nRows = matrix.rows(); - const auto nCols = matrix.cols(); - const auto rowIndices = matrix.outer(); - const auto colIndices = matrix.inner(); - auto valData = matrix.data(); - - atlas_omp_parallel_for(auto i = size_t{}; i < nRows; ++i) { - for (auto dataIdx = rowIndices[i]; dataIdx < rowIndices[i + 1]; ++dataIdx) { - const auto j = size_t(colIndices[dataIdx]); - auto&& value = valData[dataIdx]; - - if constexpr( - std::is_invocable_v) { - functor(value, i, j); - } - else if constexpr(std::is_invocable_v) { - functor(value, i, j, dataIdx); - } - else { - ATLAS_NOTIMPLEMENTED; - } +void sparseMatrixForEach(const MatrixT& matrix, const Functor& functor) { + + atlas_omp_parallel_for (auto k = 0; k < matrix.outerSize(); ++k) { + for (auto it = typename MatrixT::InnerIterator(matrix, k); it; ++it) { + functor(it.value(), it.row(), it.col()); } } } @@ -70,12 +51,12 @@ void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, sparseMatrixForEach(matrix, [&](const auto& weight, auto i, auto j) { - constexpr auto rank = std::decay_t::rank(); + 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()); mappingFunctor(weight, sourceSlice, targetSlice); - } + } else if constexpr(rank == 3) { const auto iterationFuctor = [&](auto&& sourceVars, auto&& targetVars) { mappingFunctor(weight, sourceVars, targetVars); @@ -86,7 +67,7 @@ void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, targetView.slice(i, array::Range::all(), array::Range::all()); array::helpers::ArrayForEach<0>::apply( std::tie(sourceSlice, targetSlice), iterationFuctor); - } + } else { ATLAS_NOTIMPLEMENTED; } @@ -95,14 +76,14 @@ void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, } // namespace -void ParallelTransport::do_setup(const Grid& source, const Grid& target, +void SphericalVector::do_setup(const Grid& source, const Grid& target, const Cache&) { ATLAS_NOTIMPLEMENTED; } -void ParallelTransport::do_setup(const FunctionSpace& source, +void SphericalVector::do_setup(const FunctionSpace& source, const FunctionSpace& target) { - ATLAS_TRACE("interpolation::method::ParallelTransport::do_setup"); + ATLAS_TRACE("interpolation::method::SphericalVector::do_setup"); source_ = source; target_ = target; @@ -114,20 +95,23 @@ void ParallelTransport::do_setup(const FunctionSpace& source, Interpolation(interpolationScheme_, source_, target_); setMatrix(MatrixCache(baseInterpolator)); - // Get matrix dimensions. + // Get matrix data. const auto nRows = matrix().rows(); const auto nCols = matrix().cols(); const auto nNonZeros = matrix().nonZeros(); - auto weightsReal = std::vector(nNonZeros); - auto weightsImag = std::vector(nNonZeros); + realWeights_ = + std::make_shared(nRows, nCols, nNonZeros, matrix().outer(), + matrix().inner(), matrix().data()); + + complexWeights_ = std::make_shared(nRows, nCols); + auto complexTriplets = ComplexTriplets(nNonZeros); - const auto sourceLonLats = array::make_view(source_.lonlat()); - const auto targetLonLats = array::make_view(target_.lonlat()); + sparseMatrixForEach(*realWeights_, [&](const auto& weight, auto i, auto j) { + + const auto sourceLonLats = array::make_view(source_.lonlat()); + const auto targetLonLats = array::make_view(target_.lonlat()); - // Make complex weights (would be nice if we could have a complex matrix). - sparseMatrixForEach(matrix(), - [&](auto&& weight, auto i, auto j, auto dataIdx) { const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); const auto targetLonLat = @@ -138,26 +122,23 @@ void ParallelTransport::do_setup(const FunctionSpace& source, auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - weightsReal[dataIdx] = {i, j, weight * std::cos(deltaAlpha)}; - weightsImag[dataIdx] = {i, j, weight * std::sin(deltaAlpha)}; - }); + auto idx = &weight - realWeights_->valuePtr(); - // Deal with slightly old fashioned Matrix interface - const auto buildMatrix = [&](auto& matrix, const auto& weights) { - auto tempMatrix = Matrix(nRows, nCols, weights); - matrix.swap(tempMatrix); - }; + complexTriplets[idx] = {i, j, std::polar(weight, deltaAlpha)}; + }); + complexWeights_->setFromTriplets(complexTriplets.begin(), + complexTriplets.end()); - buildMatrix(matrixReal_, weightsReal); - buildMatrix(matrixImag_, weightsImag); + ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); + ATLAS_ASSERT(realWeights_->nonZeros() == matrix().nonZeros()); } -void ParallelTransport::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } +void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } -void ParallelTransport::do_execute(const FieldSet& sourceFieldSet, +void SphericalVector::do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, Metadata& metadata) const { - ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::do_execute()"); + ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); const auto nFields = sourceFieldSet.size(); ATLAS_ASSERT(nFields == targetFieldSet.size()); @@ -167,9 +148,9 @@ void ParallelTransport::do_execute(const FieldSet& sourceFieldSet, } } -void ParallelTransport::do_execute(const Field& sourceField, Field& targetField, +void SphericalVector::do_execute(const Field& sourceField, Field& targetField, Metadata&) const { - ATLAS_TRACE("atlas::interpolation::method::ParallelTransport::do_execute()"); + ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); if (!(sourceField.variables() == 2 || sourceField.variables() == 3)) { @@ -197,7 +178,7 @@ void ParallelTransport::do_execute(const Field& sourceField, Field& targetField, } template -void ParallelTransport::interpolate_vector_field(const Field& sourceField, +void SphericalVector::interpolate_vector_field(const Field& sourceField, Field& targetField) const { if (sourceField.rank() == 2) { interpolate_vector_field(sourceField, targetField); @@ -209,32 +190,31 @@ void ParallelTransport::interpolate_vector_field(const Field& sourceField, } template -void ParallelTransport::interpolate_vector_field(const Field& sourceField, +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.); - // Matrix multiplication split in two to simulate complex variable - // multiplication. - matrixMultiply(matrixReal_, sourceView, targetView, - [](const auto& weight, auto&& sourceVars, auto&& targetVars) { - targetVars(0) += weight * sourceVars(0); - targetVars(1) += weight * sourceVars(1); - }); - matrixMultiply(matrixImag_, sourceView, targetView, - [](const auto& weight, auto&& sourceVars, auto&& targetVars) { - targetVars(0) -= weight * sourceVars(1); - targetVars(1) += weight * sourceVars(0); - }); + const auto horizontalComponent = [](const auto& weight, auto&& sourceVars, + auto&& targetVars) { + const auto targetVector = + weight * std::complex(sourceVars(0), sourceVars(1)); - if (sourceField.variables() == 3) { - matrixMultiply( - matrix(), sourceView, targetView, - [](const auto& weight, auto&& sourceVars, - auto&& targetVars) { targetVars(2) = weight * sourceVars(2); }); - } + targetVars(0) += targetVector.real(); + targetVars(1) += targetVector.imag(); + }; + + const auto verticalComponent = []( + const auto& weight, auto&& sourceVars, + auto&& targetVars) { targetVars(2) += weight * sourceVars(2); }; + + matrixMultiply(*complexWeights_, sourceView, targetView, horizontalComponent); + + if (sourceField.variables() == 2) return; + + matrixMultiply(*realWeights_, sourceView, targetView, verticalComponent); } } // namespace method diff --git a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h similarity index 73% rename from src/atlas/interpolation/method/paralleltransport/ParallelTransport.h rename to src/atlas/interpolation/method/sphericalvector/SphericalVector.h index fbc427104..a4d3703a2 100644 --- a/src/atlas/interpolation/method/paralleltransport/ParallelTransport.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -7,6 +7,11 @@ #pragma once +#include +#include + +#include + #include "atlas/functionspace/FunctionSpace.h" #include "atlas/interpolation/method/Method.h" #include "atlas/linalg/sparse.h" @@ -16,14 +21,19 @@ namespace atlas { namespace interpolation { namespace method { -class ParallelTransport : public Method { +using RealMatrix = Eigen::SparseMatrix; +using RealMatrixMap = Eigen::Map; +using ComplexMatrix = Eigen::SparseMatrix, Eigen::RowMajor>; +using ComplexTriplets = std::vector>>; + +class SphericalVector : public Method { public: - ParallelTransport(const Config& config): Method(config) { + SphericalVector(const Config& config): Method(config) { const auto& conf = dynamic_cast(config); interpolationScheme_ = conf.getSubConfiguration("scheme"); } - virtual ~ParallelTransport() override {} + virtual ~SphericalVector() override {} void print(std::ostream&) const override; const FunctionSpace& source() const override { return source_; } @@ -50,10 +60,8 @@ class ParallelTransport : public Method { FunctionSpace source_; FunctionSpace target_; - // Complex interpolation weights. We treat a (u, v) pair as a complex variable - // with u on the real number line and v on the imaginary number line. - Matrix matrixReal_; - Matrix matrixImag_; + std::shared_ptr realWeights_; + std::shared_ptr complexWeights_; }; diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index c680db212..dc4df9e29 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -81,9 +81,10 @@ ecbuild_add_test( TARGET atlas_test_interpolation_non_linear ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -ecbuild_add_test( TARGET atlas_test_interpolation_parallel_transport +ecbuild_add_test( TARGET atlas_test_interpolation_spherical_vector + CONDITION eckit_EIGEN_FOUND OMP 4 - SOURCES test_interpolation_parallel_transport.cc + SOURCES test_interpolation_spherical_vector.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) diff --git a/src/tests/interpolation/test_interpolation_parallel_transport.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc similarity index 95% rename from src/tests/interpolation/test_interpolation_parallel_transport.cc rename to src/tests/interpolation/test_interpolation_spherical_vector.cc index 23967c6f8..7e09cdf02 100644 --- a/src/tests/interpolation/test_interpolation_parallel_transport.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -154,13 +154,13 @@ CASE("cubed sphere vector interpolation") { const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear").set("adjoint", true); - const auto interpScheme = util::Config("type", "parallel-transport") + const auto interpScheme = util::Config("type", "spherical-vector") .set("scheme", baseInterpScheme); const auto cubedSphereConf = Config("source_grid", "CS-LFR-48") .set("source_mesh", "cubedsphere_dual") .set("target_grid", "O48") .set("target_mesh", "structured") - .set("file_id", "parallel_transport_cs") + .set("file_id", "spherical_vector_cs") .set("scheme", interpScheme); testInterpolation((cubedSphereConf)); @@ -170,13 +170,13 @@ CASE("finite element vector interpolation") { const auto baseInterpScheme = util::Config("type", "finite-element").set("adjoint", true); - const auto interpScheme = util::Config("type", "parallel-transport") + const auto interpScheme = util::Config("type", "spherical-vector") .set("scheme", baseInterpScheme); const auto cubedSphereConf = Config("source_grid", "O48") .set("source_mesh", "structured") .set("target_grid", "CS-LFR-48") .set("target_mesh", "cubedsphere_dual") - .set("file_id", "parallel_transport_fe") + .set("file_id", "spherical_vector_fe") .set("scheme", interpScheme); testInterpolation((cubedSphereConf)); From a6d0df97eeb8e72fb349c55d03e9979f32f44895 Mon Sep 17 00:00:00 2001 From: Oliver Lomax Date: Tue, 28 Nov 2023 19:26:51 +0000 Subject: [PATCH 06/26] Update SphericalVector.cc make_view calls erroneously in loop body. --- .../interpolation/method/sphericalvector/SphericalVector.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 4376a98ae..4a75fb09f 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -107,11 +107,11 @@ void SphericalVector::do_setup(const FunctionSpace& source, complexWeights_ = std::make_shared(nRows, nCols); auto complexTriplets = ComplexTriplets(nNonZeros); + const auto sourceLonLats = array::make_view(source_.lonlat()); + const auto targetLonLats = array::make_view(target_.lonlat()); + sparseMatrixForEach(*realWeights_, [&](const auto& weight, auto i, auto j) { - const auto sourceLonLats = array::make_view(source_.lonlat()); - const auto targetLonLats = array::make_view(target_.lonlat()); - const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); const auto targetLonLat = From 5064f449a59dd113888c07b96cce3f388c32a3b9 Mon Sep 17 00:00:00 2001 From: odlomax Date: Wed, 29 Nov 2023 15:23:36 +0000 Subject: [PATCH 07/26] Tidied up SphericalVector class. --- src/atlas/CMakeLists.txt | 2 +- .../interpolation/method/MethodFactory.cc | 6 ++- .../method/sphericalvector/SphericalVector.cc | 49 ++++++++++++------- .../method/sphericalvector/SphericalVector.h | 31 ++++++------ src/tests/interpolation/CMakeLists.txt | 19 ++++--- 5 files changed, 63 insertions(+), 44 deletions(-) diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 4a18dec2a..a98a9b680 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -669,7 +669,7 @@ interpolation/nonlinear/NonLinear.cc interpolation/nonlinear/NonLinear.h ) -if (eckit_EIGEN_FOUND) +if (atlas_HAVE_EIGEN) list (APPEND atlas_interpolation_srcs interpolation/method/sphericalvector/SphericalVector.h diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index 9f2bff494..ab078f74b 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -12,7 +12,6 @@ // for static linking #include "cubedsphere/CubedSphereBilinear.h" -#include "sphericalvector/SphericalVector.h" #include "knn/GridBoxAverage.h" #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" @@ -26,6 +25,9 @@ #include "unstructured/FiniteElement.h" #include "unstructured/UnstructuredBilinearLonLat.h" +#if ATLAS_HAVE_EIGEN +#include "sphericalvector/SphericalVector.h" +#endif namespace atlas { namespace interpolation { @@ -48,7 +50,9 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); +#if ATLAS_HAVE_EIGEN MethodBuilder(); +#endif } } link; } diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 4376a98ae..09d6b1b90 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -31,9 +31,24 @@ namespace atlas { namespace interpolation { namespace method { +using Complex = SphericalVector::Complex; + +template +using SparseMatrix = SphericalVector::SparseMatrix; +using RealMatrixMap = Eigen::Map>; +using ComplexTriplets = std::vector>; +using EckitMatrix = eckit::linalg::SparseMatrix; + namespace { + MethodBuilder __builder("spherical-vector"); +RealMatrixMap mapBaseMatrix(const EckitMatrix& baseMatrix) { + return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(), + baseMatrix.nonZeros(), baseMatrix.outer(), + baseMatrix.inner(), baseMatrix.data()); +} + template void sparseMatrixForEach(const MatrixT& matrix, const Functor& functor) { @@ -99,18 +114,15 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto nRows = matrix().rows(); const auto nCols = matrix().cols(); const auto nNonZeros = matrix().nonZeros(); - - realWeights_ = - std::make_shared(nRows, nCols, nNonZeros, matrix().outer(), - matrix().inner(), matrix().data()); + const auto realWeights = mapBaseMatrix(matrix()); complexWeights_ = std::make_shared(nRows, nCols); auto complexTriplets = ComplexTriplets(nNonZeros); - sparseMatrixForEach(*realWeights_, [&](const auto& weight, auto i, auto j) { + const auto sourceLonLats = array::make_view(source_.lonlat()); + const auto targetLonLats = array::make_view(target_.lonlat()); - const auto sourceLonLats = array::make_view(source_.lonlat()); - const auto targetLonLats = array::make_view(target_.lonlat()); + sparseMatrixForEach(realWeights, [&](const auto& weight, auto i, auto j) { const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); @@ -119,10 +131,10 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto alpha = util::greatCircleCourse(sourceLonLat, targetLonLat); - auto deltaAlpha = + const auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - auto idx = &weight - realWeights_->valuePtr(); + const auto idx = &weight - realWeights.valuePtr(); complexTriplets[idx] = {i, j, std::polar(weight, deltaAlpha)}; }); @@ -130,7 +142,6 @@ void SphericalVector::do_setup(const FunctionSpace& source, complexTriplets.end()); ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); - ATLAS_ASSERT(realWeights_->nonZeros() == matrix().nonZeros()); } void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } @@ -191,7 +202,7 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, template void SphericalVector::interpolate_vector_field(const Field& sourceField, - Field& targetField) const { + Field& targetField) const { const auto sourceView = array::make_view(sourceField); auto targetView = array::make_view(targetField); @@ -199,22 +210,22 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, const auto horizontalComponent = [](const auto& weight, auto&& sourceVars, auto&& targetVars) { - const auto targetVector = - weight * std::complex(sourceVars(0), sourceVars(1)); - + const auto sourceVector = Complex(sourceVars(0), sourceVars(1)); + const auto targetVector = weight * sourceVector; targetVars(0) += targetVector.real(); targetVars(1) += targetVector.imag(); }; - const auto verticalComponent = []( - const auto& weight, auto&& sourceVars, - auto&& targetVars) { targetVars(2) += weight * sourceVars(2); }; - matrixMultiply(*complexWeights_, sourceView, targetView, horizontalComponent); if (sourceField.variables() == 2) return; - matrixMultiply(*realWeights_, sourceView, targetView, verticalComponent); + const auto verticalComponent = []( + const auto& weight, auto&& sourceVars, + auto&& targetVars) { targetVars(2) += weight * sourceVars(2); }; + + const auto realWeights = mapBaseMatrix(matrix()); + matrixMultiply(realWeights, sourceView, targetView, verticalComponent); } } // namespace method diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index a4d3703a2..1e4a9a715 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -21,17 +21,18 @@ namespace atlas { namespace interpolation { namespace method { -using RealMatrix = Eigen::SparseMatrix; -using RealMatrixMap = Eigen::Map; -using ComplexMatrix = Eigen::SparseMatrix, Eigen::RowMajor>; -using ComplexTriplets = std::vector>>; - class SphericalVector : public Method { public: - SphericalVector(const Config& config): Method(config) { + + using Complex = std::complex; + + template + using SparseMatrix = Eigen::SparseMatrix; + using ComplexMatrix = SparseMatrix; + + SphericalVector(const Config& config) : Method(config) { const auto& conf = dynamic_cast(config); interpolationScheme_ = conf.getSubConfiguration("scheme"); - } virtual ~SphericalVector() override {} @@ -39,20 +40,22 @@ class SphericalVector : public Method { 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(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, + Metadata& metadata) const override; + void do_execute(const Field& sourceField, Field& targetField, + Metadata& metadata) const override; private: - template + template void interpolate_vector_field(const Field& source, Field& target) const; - template + 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 FunctionSpace& source, + const FunctionSpace& target) override; void do_setup(const Grid& source, const Grid& target, const Cache&) override; eckit::LocalConfiguration interpolationScheme_; @@ -60,9 +63,7 @@ class SphericalVector : public Method { FunctionSpace source_; FunctionSpace target_; - std::shared_ptr realWeights_; std::shared_ptr complexWeights_; - }; } // namespace method diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index dc4df9e29..429323cd4 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -81,14 +81,6 @@ ecbuild_add_test( TARGET atlas_test_interpolation_non_linear ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -ecbuild_add_test( TARGET atlas_test_interpolation_spherical_vector - CONDITION eckit_EIGEN_FOUND - OMP 4 - SOURCES test_interpolation_spherical_vector.cc - LIBS atlas - ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} -) - ecbuild_add_test( TARGET atlas_test_interpolation_bilinear COMMAND atlas_test_interpolation_structured2D ARGS --scheme linear ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} @@ -130,4 +122,15 @@ ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +if ( atlas_HAVE_EIGEN ) + +ecbuild_add_test( TARGET atlas_test_interpolation_spherical_vector + OMP 4 + SOURCES test_interpolation_spherical_vector.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + +endif() + endif() From 6b43f8452a95d52f83946d6bb4b705d7f3d4cedf Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 30 Nov 2023 11:09:23 +0000 Subject: [PATCH 08/26] Minor cosmetic changes. --- .../interpolation/method/sphericalvector/SphericalVector.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 2c5ce55ad..7d2ad7d2b 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -134,7 +134,7 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - const auto idx = &weight - realWeights.valuePtr(); + const auto idx = std::distance(realWeights.valuePtr(), &weight); complexTriplets[idx] = {int(i), int(j), std::polar(weight, deltaAlpha)}; }); @@ -163,7 +163,7 @@ void SphericalVector::do_execute(const Field& sourceField, Field& targetField, Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); - const auto fieldType = targetField.metadata().getString("type", ""); + const auto fieldType = sourceField.metadata().getString("type", ""); if (fieldType != "vector") { auto metadata = Metadata(); From 35611a8726637a73a3e2b13670c691f5f87f16fd Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 30 Nov 2023 11:33:19 +0000 Subject: [PATCH 09/26] Replaced optional compilation with #if ATLAS_HAVE_EIGEN in source file and header. --- src/atlas/CMakeLists.txt | 11 ++--------- src/atlas/interpolation/method/MethodFactory.cc | 8 ++++---- .../method/sphericalvector/SphericalVector.cc | 5 +++++ .../method/sphericalvector/SphericalVector.h | 5 +++++ 4 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index a98a9b680..4448e3fd0 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -632,6 +632,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 @@ -669,15 +671,6 @@ interpolation/nonlinear/NonLinear.cc interpolation/nonlinear/NonLinear.h ) -if (atlas_HAVE_EIGEN) - -list (APPEND atlas_interpolation_srcs -interpolation/method/sphericalvector/SphericalVector.h -interpolation/method/sphericalvector/SphericalVector.cc -) - -endif() - list( APPEND atlas_linalg_srcs linalg/Indexing.h linalg/Introspection.h diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index ab078f74b..cc30db2bd 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -8,6 +8,7 @@ * nor does it submit to any jurisdiction. */ +#include "atlas/library/defines.h" #include "MethodFactory.h" // for static linking @@ -16,6 +17,9 @@ #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" #include "knn/NearestNeighbour.h" +#if ATLAS_HAVE_EIGEN +#include "sphericalvector/SphericalVector.h" +#endif #include "structured/Cubic2D.h" #include "structured/Cubic3D.h" #include "structured/Linear2D.h" @@ -25,10 +29,6 @@ #include "unstructured/FiniteElement.h" #include "unstructured/UnstructuredBilinearLonLat.h" -#if ATLAS_HAVE_EIGEN -#include "sphericalvector/SphericalVector.h" -#endif - namespace atlas { namespace interpolation { diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 7d2ad7d2b..81450bb33 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -5,6 +5,9 @@ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ +#include "atlas/library/defines.h" +#if ATLAS_HAVE_EIGEN + #include #include @@ -237,3 +240,5 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, } // namespace method } // namespace interpolation } // namespace atlas + +#endif diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index de9b65f95..780c43a79 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -5,6 +5,9 @@ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ +#include "atlas/library/defines.h" +#if ATLAS_HAVE_EIGEN + #pragma once #include @@ -92,3 +95,5 @@ class SphericalVector : public Method { } // namespace method } // namespace interpolation } // namespace atlas + +#endif From a913b180f744c8e8b9107080860f3e934582667f Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 30 Nov 2023 11:38:41 +0000 Subject: [PATCH 10/26] Removed redundant macros. --- src/atlas/interpolation/method/MethodFactory.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index cc30db2bd..f68ace367 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -17,9 +17,7 @@ #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" #include "knn/NearestNeighbour.h" -#if ATLAS_HAVE_EIGEN #include "sphericalvector/SphericalVector.h" -#endif #include "structured/Cubic2D.h" #include "structured/Cubic3D.h" #include "structured/Linear2D.h" From c4481aa1e9fc1305c9c2e62a0a5259057fc4dbf7 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 30 Nov 2023 15:07:11 +0000 Subject: [PATCH 11/26] Removed static factory linking. --- src/atlas/interpolation/method/MethodFactory.cc | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index f68ace367..de1739863 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -8,7 +8,6 @@ * nor does it submit to any jurisdiction. */ -#include "atlas/library/defines.h" #include "MethodFactory.h" // for static linking @@ -17,7 +16,6 @@ #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" @@ -48,9 +46,6 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); -#if ATLAS_HAVE_EIGEN - MethodBuilder(); -#endif } } link; } From 6107e913e6274985bc2753999da75162b00b8504 Mon Sep 17 00:00:00 2001 From: odlomax Date: Fri, 1 Dec 2023 18:28:57 +0000 Subject: [PATCH 12/26] Fused horizontal and vertical component matrix-multiplications. TODO: Write test for 3-vector field (note: this works when hacking current 2-vector test). --- .../method/sphericalvector/SphericalVector.cc | 71 +++++++++++-------- .../test_interpolation_spherical_vector.cc | 10 +-- 2 files changed, 47 insertions(+), 34 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 81450bb33..e5ac80fb2 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -10,6 +10,7 @@ #include #include +#include #include "atlas/array/ArrayView.h" #include "atlas/array/helpers/ArrayForEach.h" @@ -55,9 +56,9 @@ RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { template void sparseMatrixForEach(const MatrixT& matrix, const Functor& functor) { - atlas_omp_parallel_for (auto k = 0; k < matrix.outerSize(); ++k) { + atlas_omp_parallel_for(auto k = 0; k < matrix.outerSize(); ++k) { for (auto it = typename MatrixT::InnerIterator(matrix, k); it; ++it) { - functor(it.value(), it.row(), it.col()); + functor(it.value(), it.row(), it.col()); } } } @@ -65,27 +66,25 @@ void sparseMatrixForEach(const MatrixT& matrix, const Functor& functor) { template void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, - TargetView&& targetView, const Functor& mappingFunctor) { + TargetView&& targetView, const Functor& multiplyFunctor) { sparseMatrixForEach(matrix, [&](const auto& weight, auto i, auto j) { - constexpr auto rank = std::decay_t::rank(); - if constexpr(rank == 2) { + 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()); - mappingFunctor(weight, sourceSlice, targetSlice); - } - else if constexpr(rank == 3) { - const auto iterationFuctor = [&](auto&& sourceVars, auto&& targetVars) { - mappingFunctor(weight, sourceVars, targetVars); - }; + multiplyFunctor(weight, sourceSlice, targetSlice); + } + 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), iterationFuctor); - } + std::tie(sourceSlice, targetSlice), + [&](auto&&... slices) { multiplyFunctor(weight, slices...); }); + } else { ATLAS_NOTIMPLEMENTED; } @@ -95,12 +94,12 @@ void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, } // namespace void SphericalVector::do_setup(const Grid& source, const Grid& target, - const Cache&) { + const Cache&) { ATLAS_NOTIMPLEMENTED; } void SphericalVector::do_setup(const FunctionSpace& source, - const FunctionSpace& target) { + const FunctionSpace& target) { ATLAS_TRACE("interpolation::method::SphericalVector::do_setup"); source_ = source; target_ = target; @@ -109,9 +108,7 @@ void SphericalVector::do_setup(const FunctionSpace& source, return; } - const auto baseInterpolator = - Interpolation(interpolationScheme_, source_, target_); - setMatrix(MatrixCache(baseInterpolator)); + setMatrix(Interpolation(interpolationScheme_, source_, target_)); // Get matrix data. const auto nRows = matrix().rows(); @@ -150,8 +147,8 @@ void SphericalVector::do_setup(const FunctionSpace& source, void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } void SphericalVector::do_execute(const FieldSet& sourceFieldSet, - FieldSet& targetFieldSet, - Metadata& metadata) const { + FieldSet& targetFieldSet, + Metadata& metadata) const { ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); const auto nFields = sourceFieldSet.size(); @@ -163,7 +160,7 @@ void SphericalVector::do_execute(const FieldSet& sourceFieldSet, } void SphericalVector::do_execute(const Field& sourceField, Field& targetField, - Metadata&) const { + Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); const auto fieldType = sourceField.metadata().getString("type", ""); @@ -199,7 +196,7 @@ void SphericalVector::do_execute(const Field& sourceField, Field& targetField, template void SphericalVector::interpolate_vector_field(const Field& sourceField, - Field& targetField) const { + Field& targetField) const { if (sourceField.rank() == 2) { interpolate_vector_field(sourceField, targetField); } else if (sourceField.rank() == 3) { @@ -225,16 +222,32 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, targetVars(1) += targetVector.imag(); }; - matrixMultiply(*complexWeights_, sourceView, targetView, horizontalComponent); + if (sourceField.variables() == 2) { + matrixMultiply(*complexWeights_, sourceView, targetView, + horizontalComponent); + return; + } else if (sourceField.variables() == 3) { - if (sourceField.variables() == 2) return; + const auto magnitudesArray = matrix().data(); + const auto* weightsBegin = complexWeights_->valuePtr(); + const auto weightMagnitude = [&](const auto& weight) { + const auto idx = std::distance(weightsBegin, &weight); + return magnitudesArray[idx]; + }; - const auto verticalComponent = []( - const auto& weight, auto&& sourceVars, - auto&& targetVars) { targetVars(2) += weight * sourceVars(2); }; + const auto horizontalAndVerticalComponent = [&]( + const auto& weight, auto&& sourceVars, auto&& targetVars) { + horizontalComponent(weight, sourceVars, targetVars); + targetVars(2) += weightMagnitude(weight) * sourceVars(2); + }; - const auto realWeights = makeMatrixMap(matrix()); - matrixMultiply(realWeights, sourceView, targetView, verticalComponent); + matrixMultiply(*complexWeights_, sourceView, targetView, + horizontalAndVerticalComponent); + + return; + } + + ATLAS_NOTIMPLEMENTED; } } // namespace method diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index 81a5618e3..7148616d6 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -86,19 +86,18 @@ void testMeshInterpolation(const Config& config) { auto sourceField = array::make_view( sourceFieldSet.add(sourceFunctionSpace.createField( option::name("test field") | option::levels(1) | - option::variables(3) | option::type("vector")))); + option::variables(2) | option::type("vector")))); auto targetField = array::make_view( targetFieldSet.add(targetFunctionSpace.createField( option::name("test field") | option::levels(1) | - option::variables(3) | option::type("vector")))); + option::variables(2) | option::type("vector")))); ArrayForEach<0>::apply(std::tie(sourceLonLat, sourceField), [](auto&& lonLat, auto&& sourceColumn) { ArrayForEach<0>::apply(std::tie(sourceColumn), [&](auto&& sourceElem) { std::tie(sourceElem(0), sourceElem(1)) = vortexField(lonLat(0), lonLat(1)); - sourceElem(2) = 0.; }); }); @@ -137,7 +136,7 @@ void testMeshInterpolation(const Config& config) { CASE("cubed sphere vector interpolation") { const auto baseInterpScheme = - util::Config("type", "cubedsphere-bilinear").set("adjoint", true); + util::Config("type", "cubedsphere-bilinear"); const auto interpScheme = util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); const auto cubedSphereConf = Config("source_grid", "CS-LFR-48") @@ -154,7 +153,7 @@ CASE("cubed sphere vector interpolation") { CASE("finite element vector interpolation") { const auto baseInterpScheme = - util::Config("type", "finite-element").set("adjoint", true); + util::Config("type", "finite-element"); const auto interpScheme = util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); const auto cubedSphereConf = Config("source_grid", "O48") @@ -167,6 +166,7 @@ CASE("finite element vector interpolation") { testMeshInterpolation((cubedSphereConf)); } + } } From 77856f3f915c0f57927faf8efe63948e1b19ae24 Mon Sep 17 00:00:00 2001 From: odlomax Date: Fri, 1 Dec 2023 23:51:57 +0000 Subject: [PATCH 13/26] Tidied fused loop. --- .../method/sphericalvector/SphericalVector.cc | 69 +++++++++++-------- 1 file changed, 40 insertions(+), 29 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index e5ac80fb2..1099877e9 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -58,34 +58,49 @@ void sparseMatrixForEach(const MatrixT& matrix, const Functor& functor) { atlas_omp_parallel_for(auto k = 0; k < matrix.outerSize(); ++k) { for (auto it = typename MatrixT::InnerIterator(matrix, k); it; ++it) { - functor(it.value(), it.row(), it.col()); + functor(it.row(), it.col(), it.value()); } } } -template -void matrixMultiply(const MatrixT& matrix, SourceView&& sourceView, - TargetView&& targetView, const Functor& multiplyFunctor) { +template +void sparseMatrixForEach(const MatrixT1& matrix1, const MatrixT2& matrix2, + const Functor& functor) { - sparseMatrixForEach(matrix, [&](const auto& weight, auto i, auto j) { + atlas_omp_parallel_for(auto k = 0; k < matrix1.outerSize(); ++k) { + for (auto[it1, it2] = + std::make_pair(typename MatrixT1::InnerIterator(matrix1, k), + typename MatrixT2::InnerIterator(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) { + + sparseMatrixForEach(matrices..., [&](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(weight, sourceSlice, targetSlice); - } - else if constexpr (Rank == 3) { + 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&&... slices) { multiplyFunctor(weight, slices...); }); - } - else { + [&](auto&& sourceVars, auto&& targetVars) { + multiplyFunctor(sourceVars, targetVars, weights...); + }); + } else { ATLAS_NOTIMPLEMENTED; } }); @@ -122,7 +137,7 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto sourceLonLats = array::make_view(source_.lonlat()); const auto targetLonLats = array::make_view(target_.lonlat()); - sparseMatrixForEach(realWeights, [&](const auto& weight, auto i, auto j) { + sparseMatrixForEach(realWeights, [&](auto i, auto j, const auto& weight) { const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); @@ -214,35 +229,31 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, auto targetView = array::make_view(targetField); targetView.assign(0.); - const auto horizontalComponent = [](const auto& weight, auto&& sourceVars, - auto&& targetVars) { + const auto horizontalComponent = [](const auto& sourceVars, auto& targetVars, + const auto& complexWeight) { const auto sourceVector = Complex(sourceVars(0), sourceVars(1)); - const auto targetVector = weight * sourceVector; + const auto targetVector = complexWeight * sourceVector; targetVars(0) += targetVector.real(); targetVars(1) += targetVector.imag(); }; if (sourceField.variables() == 2) { - matrixMultiply(*complexWeights_, sourceView, targetView, - horizontalComponent); + matrixMultiply(sourceView, targetView, horizontalComponent, + *complexWeights_); return; } else if (sourceField.variables() == 3) { - const auto magnitudesArray = matrix().data(); - const auto* weightsBegin = complexWeights_->valuePtr(); - const auto weightMagnitude = [&](const auto& weight) { - const auto idx = std::distance(weightsBegin, &weight); - return magnitudesArray[idx]; - }; + const auto realWeights = makeMatrixMap(matrix()); const auto horizontalAndVerticalComponent = [&]( - const auto& weight, auto&& sourceVars, auto&& targetVars) { - horizontalComponent(weight, sourceVars, targetVars); - targetVars(2) += weightMagnitude(weight) * sourceVars(2); + const auto& sourceVars, auto& targetVars, const auto& complexWeight, + const auto& realWeight) { + horizontalComponent(sourceVars, targetVars, complexWeight); + targetVars(2) += realWeight * sourceVars(2); }; - matrixMultiply(*complexWeights_, sourceView, targetView, - horizontalAndVerticalComponent); + matrixMultiply(sourceView, targetView, horizontalAndVerticalComponent, + *complexWeights_, realWeights); return; } From 346f7f1c904abfb0bd38590c249a565ed4177563 Mon Sep 17 00:00:00 2001 From: odlomax Date: Tue, 5 Dec 2023 10:35:14 +0000 Subject: [PATCH 14/26] Uncovered and fixed differences in eckit and Eigen3 CRS format. Also added tests for 3-vector and 2 vector fields. --- .../method/sphericalvector/SphericalVector.cc | 54 +++-- .../method/sphericalvector/SphericalVector.h | 2 + .../test_interpolation_cubedsphere.cc | 6 +- .../test_interpolation_spherical_vector.cc | 198 ++++++++++++------ 4 files changed, 170 insertions(+), 90 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 1099877e9..4f8e643d3 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -41,6 +41,7 @@ template using SparseMatrix = SphericalVector::SparseMatrix; using RealMatrixMap = Eigen::Map>; using ComplexTriplets = std::vector>; +using RealTriplets = std::vector>; using EckitMatrix = eckit::linalg::SparseMatrix; namespace { @@ -53,24 +54,27 @@ RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { baseMatrix.inner(), baseMatrix.data()); } -template -void sparseMatrixForEach(const MatrixT& matrix, const Functor& functor) { +template +auto getInnerIt(const Matrix& matrix, int k) { + return typename Matrix::InnerIterator(matrix, k); +} + +template +void sparseMatrixForEach(const Functor& functor, const Matrix& matrix) { atlas_omp_parallel_for(auto k = 0; k < matrix.outerSize(); ++k) { - for (auto it = typename MatrixT::InnerIterator(matrix, k); it; ++it) { + for (auto it = getInnerIt(matrix, k); it; ++it) { functor(it.row(), it.col(), it.value()); } } } -template -void sparseMatrixForEach(const MatrixT1& matrix1, const MatrixT2& matrix2, - const Functor& functor) { - +template +void sparseMatrixForEach(const Functor& functor, const Matrix1& matrix1, + const Matrix2& matrix2) { atlas_omp_parallel_for(auto k = 0; k < matrix1.outerSize(); ++k) { - for (auto[it1, it2] = - std::make_pair(typename MatrixT1::InnerIterator(matrix1, k), - typename MatrixT2::InnerIterator(matrix2, 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()); } @@ -83,14 +87,13 @@ void matrixMultiply(const SourceView& sourceView, TargetView& targetView, const Functor& multiplyFunctor, const Matrices&... matrices) { - sparseMatrixForEach(matrices..., [&](auto i, auto j, const auto&... weights) { - + 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) { + } else if constexpr(Rank == 3) { const auto sourceSlice = sourceView.slice(j, array::Range::all(), array::Range::all()); auto targetSlice = @@ -103,7 +106,9 @@ void matrixMultiply(const SourceView& sourceView, TargetView& targetView, } else { ATLAS_NOTIMPLEMENTED; } - }); + }; + + sparseMatrixForEach(multiplyColumn, matrices...); } } // namespace @@ -129,16 +134,19 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto nRows = matrix().rows(); const auto nCols = matrix().cols(); const auto nNonZeros = matrix().nonZeros(); - const auto realWeights = makeMatrixMap(matrix()); + 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()); - sparseMatrixForEach(realWeights, [&](auto i, auto j, const auto& weight) { - + const auto setComplexWeights = [&](auto i, auto j, const auto& weight) { const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); const auto targetLonLat = @@ -149,12 +157,16 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - const auto idx = std::distance(realWeights.valuePtr(), &weight); + const auto idx = std::distance(baseWeights.valuePtr(), &weight); complexTriplets[idx] = {int(i), int(j), std::polar(weight, deltaAlpha)}; - }); + realTriplets[idx] = {int(i), int(j), weight}; + }; + + sparseMatrixForEach(setComplexWeights, baseWeights); complexWeights_->setFromTriplets(complexTriplets.begin(), complexTriplets.end()); + realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end()); ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); } @@ -243,8 +255,6 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, return; } else if (sourceField.variables() == 3) { - const auto realWeights = makeMatrixMap(matrix()); - const auto horizontalAndVerticalComponent = [&]( const auto& sourceVars, auto& targetVars, const auto& complexWeight, const auto& realWeight) { @@ -253,7 +263,7 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, }; matrixMultiply(sourceView, targetView, horizontalAndVerticalComponent, - *complexWeights_, realWeights); + *complexWeights_, *realWeights_); return; } diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index 780c43a79..f3bc5e1ac 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -31,6 +31,7 @@ class SphericalVector : public Method { template using SparseMatrix = Eigen::SparseMatrix; using ComplexMatrix = SparseMatrix; + using RealMatrix = SparseMatrix; /// @brief Interpolation post-processor for vector field interpolation /// @@ -90,6 +91,7 @@ class SphericalVector : public Method { FunctionSpace target_; std::shared_ptr complexWeights_; + std::shared_ptr realWeights_; }; } // namespace method 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 index 7148616d6..143780b2b 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -27,9 +27,12 @@ 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 vortexField(double lon, double lat) { +std::pair vortexHorizontal(double lon, double lat) { // set hLon and hLat step size. const double hLon = 0.0001; @@ -38,18 +41,21 @@ std::pair vortexField(double lon, double lat) { // Get finite differences. // Set u. - const double u = (util::function::vortex_rollup(lon, lat + 0.5 * hLat, 0.1) - - util::function::vortex_rollup(lon, lat - 0.5 * hLat, 0.1)) / + 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 = - -(util::function::vortex_rollup(lon + 0.5 * hLon, lat, 0.1) - - util::function::vortex_rollup(lon - 0.5 * hLon, lat, 0.1)) / - (hLon * std::cos(lat * util::Constants::degreesToRadians())); + 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(); @@ -63,17 +69,45 @@ void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { gmsh.write(fieldSet, functionSpace); } -void testMeshInterpolation(const Config& config) { - - const auto sourceGrid = Grid(config.getString("source_grid")); - const auto sourceMesh = - MeshGenerator(config.getString("source_mesh")).generate(sourceGrid); - const auto sourceFunctionSpace = functionspace::NodeColumns(sourceMesh); +// 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); +} - const auto targetGrid = Grid(config.getString("target_grid")); - const auto targetMesh = - MeshGenerator(config.getString("target_mesh")).generate(targetGrid); - const auto targetFunctionSpace = functionspace::NodeColumns(targetMesh); +// 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{}; @@ -83,23 +117,32 @@ void testMeshInterpolation(const Config& config) { const auto targetLonLat = array::make_view(targetFunctionSpace.lonlat()); - auto sourceField = array::make_view( - sourceFieldSet.add(sourceFunctionSpace.createField( - option::name("test field") | option::levels(1) | - option::variables(2) | option::type("vector")))); + auto fieldSpec = + FieldSpecFixtures::get(config.getString("field_spec_fixture")); + if constexpr (Rank == 3) fieldSpec.set("levels", 1); - auto targetField = array::make_view( - targetFieldSet.add(targetFunctionSpace.createField( - option::name("test field") | option::levels(1) | - option::variables(2) | option::type("vector")))); + 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) { - ArrayForEach<0>::apply(std::tie(sourceColumn), [&](auto&& sourceElem) { + + const auto setElems = [&](auto&& sourceElem) { std::tie(sourceElem(0), sourceElem(1)) = - vortexField(lonLat(0), lonLat(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); @@ -107,66 +150,91 @@ void testMeshInterpolation(const Config& config) { interp.execute(sourceFieldSet, targetFieldSet); targetFieldSet.haloExchange(); - auto errorField = array::make_view( - targetFieldSet.add(targetFunctionSpace.createField( - option::name("error field") | option::levels(1)))); + 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) { - ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), - [&](auto&& targetElem, auto&& errorElem) { - - auto deltaVal = vortexField(lonLat(0), lonLat(1)); - deltaVal.first -= targetElem(0); - deltaVal.second -= targetElem(1); - - errorElem = std::sqrt(deltaVal.first * deltaVal.first + - deltaVal.second * deltaVal.second); - maxError = std::max(maxError, static_cast(errorElem)); - }); - }); + 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") { +CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { - const auto baseInterpScheme = - util::Config("type", "cubedsphere-bilinear"); + const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear"); const auto interpScheme = util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); - const auto cubedSphereConf = Config("source_grid", "CS-LFR-48") - .set("source_mesh", "cubedsphere_dual") - .set("target_grid", "O48") - .set("target_mesh", "structured") - .set("file_id", "spherical_vector_cs") + 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); - testMeshInterpolation((cubedSphereConf)); + testInterpolation((cubedSphereConf)); } -CASE("finite element vector interpolation") { +CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { - const auto baseInterpScheme = - util::Config("type", "finite-element"); + const auto baseInterpScheme = util::Config("type", "cubedsphere-bilinear"); const auto interpScheme = util::Config("type", "spherical-vector").set("scheme", baseInterpScheme); - const auto cubedSphereConf = Config("source_grid", "O48") - .set("source_mesh", "structured") - .set("target_grid", "CS-LFR-48") - .set("target_mesh", "cubedsphere_dual") - .set("file_id", "spherical_vector_fe") + 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.00015); + .set("tol", 0.00096); - testMeshInterpolation((cubedSphereConf)); + 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)); +} } } From ee4c3adb5ae998a53c11dbc50ef6d0a6ee6567e9 Mon Sep 17 00:00:00 2001 From: odlomax Date: Tue, 5 Dec 2023 10:42:08 +0000 Subject: [PATCH 15/26] Added multiple levels to 3d fields. --- src/tests/interpolation/test_interpolation_spherical_vector.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index 143780b2b..1e68b8e91 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -119,7 +119,7 @@ void testInterpolation(const Config& config) { auto fieldSpec = FieldSpecFixtures::get(config.getString("field_spec_fixture")); - if constexpr (Rank == 3) fieldSpec.set("levels", 1); + if constexpr (Rank == 3) fieldSpec.set("levels", 2); auto sourceField = array::make_view( sourceFieldSet.add(sourceFunctionSpace.createField(fieldSpec))); From a96057eee9cd3d5adc462219c90dc5250248e0da Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 7 Dec 2023 15:52:48 +0000 Subject: [PATCH 16/26] Add SphericalVector to MethodFactory --- .../interpolation/method/MethodFactory.cc | 2 ++ .../method/sphericalvector/SphericalVector.cc | 18 +++++++++++++++--- .../method/sphericalvector/SphericalVector.h | 9 ++++++--- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index de1739863..e666e2bf2 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" @@ -46,6 +47,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 index 4f8e643d3..96318f250 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -6,7 +6,6 @@ */ #include "atlas/library/defines.h" -#if ATLAS_HAVE_EIGEN #include #include @@ -37,17 +36,21 @@ namespace method { using Complex = SphericalVector::Complex; +#if ATLAS_HAVE_EIGEN template using SparseMatrix = SphericalVector::SparseMatrix; using RealMatrixMap = Eigen::Map>; using ComplexTriplets = std::vector>; using RealTriplets = std::vector>; +#endif + using EckitMatrix = eckit::linalg::SparseMatrix; namespace { MethodBuilder __builder("spherical-vector"); +#if ATLAS_HAVE_EIGEN RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(), baseMatrix.nonZeros(), baseMatrix.outer(), @@ -110,6 +113,7 @@ void matrixMultiply(const SourceView& sourceView, TargetView& targetView, sparseMatrixForEach(multiplyColumn, matrices...); } +#endif } // namespace @@ -128,6 +132,8 @@ void SphericalVector::do_setup(const FunctionSpace& source, return; } +#if ATLAS_HAVE_EIGEN + setMatrix(Interpolation(interpolationScheme_, source_, target_)); // Get matrix data. @@ -169,6 +175,10 @@ void SphericalVector::do_setup(const FunctionSpace& source, realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end()); ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); + +#else + ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen"); +#endif } void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } @@ -241,6 +251,7 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, auto targetView = array::make_view(targetField); targetView.assign(0.); +#if ATLAS_HAVE_EIGEN const auto horizontalComponent = [](const auto& sourceVars, auto& targetVars, const auto& complexWeight) { const auto sourceVector = Complex(sourceVars(0), sourceVars(1)); @@ -267,6 +278,9 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, return; } +#else + ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen"); +#endif ATLAS_NOTIMPLEMENTED; } @@ -274,5 +288,3 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, } // namespace method } // namespace interpolation } // namespace atlas - -#endif diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index f3bc5e1ac..8f22a8c62 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -6,14 +6,15 @@ */ #include "atlas/library/defines.h" -#if ATLAS_HAVE_EIGEN #pragma once #include #include +#if ATLAS_HAVE_EIGEN #include +#endif #include "atlas/functionspace/FunctionSpace.h" #include "atlas/interpolation/method/Method.h" @@ -28,10 +29,12 @@ class SphericalVector : public Method { public: using Complex = std::complex; +#if ATLAS_HAVE_EIGEN template using SparseMatrix = Eigen::SparseMatrix; using ComplexMatrix = SparseMatrix; using RealMatrix = SparseMatrix; +#endif /// @brief Interpolation post-processor for vector field interpolation /// @@ -90,12 +93,12 @@ class SphericalVector : public Method { FunctionSpace source_; FunctionSpace target_; +#if ATLAS_HAVE_EIGEN std::shared_ptr complexWeights_; std::shared_ptr realWeights_; +#endif }; } // namespace method } // namespace interpolation } // namespace atlas - -#endif From 8b3e8f4cef0d26b5667c574d9a76476779746855 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 7 Dec 2023 18:24:09 +0000 Subject: [PATCH 17/26] Added more consistent types to iteration indices. --- .../interpolation/method/sphericalvector/SphericalVector.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 96318f250..1aeab59d9 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -65,7 +65,8 @@ auto getInnerIt(const Matrix& matrix, int k) { template void sparseMatrixForEach(const Functor& functor, const Matrix& matrix) { - atlas_omp_parallel_for(auto k = 0; k < matrix.outerSize(); ++k) { + using Index = decltype (matrix.outerSize()); + 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()); } @@ -75,7 +76,8 @@ void sparseMatrixForEach(const Functor& functor, const Matrix& matrix) { template void sparseMatrixForEach(const Functor& functor, const Matrix1& matrix1, const Matrix2& matrix2) { - atlas_omp_parallel_for(auto k = 0; k < matrix1.outerSize(); ++k) { + using Index = decltype (matrix1.outerSize()); + 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) { From d63503fd54a8a1ad179395fcccb8c3d0ccb7ceab Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 7 Dec 2023 18:31:04 +0000 Subject: [PATCH 18/26] Further index consistency added. --- .../interpolation/method/sphericalvector/SphericalVector.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 1aeab59d9..b4692d791 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -57,8 +57,8 @@ RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { baseMatrix.inner(), baseMatrix.data()); } -template -auto getInnerIt(const Matrix& matrix, int k) { +template +auto getInnerIt(const Matrix& matrix, Index k) { return typename Matrix::InnerIterator(matrix, k); } From 636c6d80c1a88eb76fc5a5791d15594bbf08de09 Mon Sep 17 00:00:00 2001 From: odlomax Date: Fri, 8 Dec 2023 09:59:06 +0000 Subject: [PATCH 19/26] Removed superfluous templates. --- .../method/sphericalvector/SphericalVector.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index b4692d791..16619a497 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -57,15 +57,15 @@ RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { baseMatrix.inner(), baseMatrix.data()); } -template -auto getInnerIt(const Matrix& matrix, Index k) { +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 = decltype (matrix.outerSize()); + 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()); @@ -76,7 +76,7 @@ void sparseMatrixForEach(const Functor& functor, const Matrix& matrix) { template void sparseMatrixForEach(const Functor& functor, const Matrix1& matrix1, const Matrix2& matrix2) { - using Index = decltype (matrix1.outerSize()); + 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)); From 3d16b6b0c75d0112852dc9df9c1358ab7e78189e Mon Sep 17 00:00:00 2001 From: odlomax Date: Fri, 8 Dec 2023 16:07:23 +0000 Subject: [PATCH 20/26] Tided up macros. --- .../method/sphericalvector/SphericalVector.cc | 36 ++++++++----------- .../method/sphericalvector/SphericalVector.h | 33 ++++++++++++++--- 2 files changed, 42 insertions(+), 27 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 16619a497..4f7da3481 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -34,23 +34,24 @@ namespace atlas { namespace interpolation { namespace method { -using Complex = SphericalVector::Complex; +namespace { +MethodBuilder __builder("spherical-vector"); +} #if ATLAS_HAVE_EIGEN + +using Complex = SphericalVector::Complex; + template using SparseMatrix = SphericalVector::SparseMatrix; using RealMatrixMap = Eigen::Map>; using ComplexTriplets = std::vector>; using RealTriplets = std::vector>; -#endif using EckitMatrix = eckit::linalg::SparseMatrix; namespace { -MethodBuilder __builder("spherical-vector"); - -#if ATLAS_HAVE_EIGEN RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(), baseMatrix.nonZeros(), baseMatrix.outer(), @@ -64,7 +65,6 @@ auto getInnerIt(const Matrix& matrix, typename Matrix::Index 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) { @@ -115,7 +115,6 @@ void matrixMultiply(const SourceView& sourceView, TargetView& targetView, sparseMatrixForEach(multiplyColumn, matrices...); } -#endif } // namespace @@ -134,8 +133,6 @@ void SphericalVector::do_setup(const FunctionSpace& source, return; } -#if ATLAS_HAVE_EIGEN - setMatrix(Interpolation(interpolationScheme_, source_, target_)); // Get matrix data. @@ -154,7 +151,7 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto sourceLonLats = array::make_view(source_.lonlat()); const auto targetLonLats = array::make_view(target_.lonlat()); - const auto setComplexWeights = [&](auto i, auto j, const auto& weight) { + const auto setWeights = [&](auto i, auto j, const auto& baseWeight) { const auto sourceLonLat = PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); const auto targetLonLat = @@ -165,22 +162,19 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - const auto idx = std::distance(baseWeights.valuePtr(), &weight); + const auto idx = std::distance(baseWeights.valuePtr(), &baseWeight); - complexTriplets[idx] = {int(i), int(j), std::polar(weight, deltaAlpha)}; - realTriplets[idx] = {int(i), int(j), weight}; + complexTriplets[idx] = {int(i), int(j), std::polar(baseWeight, deltaAlpha)}; + realTriplets[idx] = {int(i), int(j), baseWeight}; }; - sparseMatrixForEach(setComplexWeights, baseWeights); + sparseMatrixForEach(setWeights, baseWeights); complexWeights_->setFromTriplets(complexTriplets.begin(), complexTriplets.end()); realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end()); ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros()); - -#else - ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen"); -#endif + ATLAS_ASSERT(realWeights_->nonZeros() == matrix().nonZeros()); } void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } @@ -253,7 +247,6 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, auto targetView = array::make_view(targetField); targetView.assign(0.); -#if ATLAS_HAVE_EIGEN const auto horizontalComponent = [](const auto& sourceVars, auto& targetVars, const auto& complexWeight) { const auto sourceVector = Complex(sourceVars(0), sourceVars(1)); @@ -280,13 +273,12 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, return; } -#else - ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen"); -#endif 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 8f22a8c62..5807fd136 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -25,16 +25,16 @@ namespace atlas { namespace interpolation { namespace method { +#if ATLAS_HAVE_EIGEN class SphericalVector : public Method { public: using Complex = std::complex; -#if ATLAS_HAVE_EIGEN template using SparseMatrix = Eigen::SparseMatrix; using ComplexMatrix = SparseMatrix; using RealMatrix = SparseMatrix; -#endif + /// @brief Interpolation post-processor for vector field interpolation /// @@ -55,7 +55,7 @@ class SphericalVector : public Method { const auto& conf = dynamic_cast(config); interpolationScheme_ = conf.getSubConfiguration("scheme"); } - virtual ~SphericalVector() override {} + ~SphericalVector() override {} void print(std::ostream&) const override; const FunctionSpace& source() const override { return source_; } @@ -93,11 +93,34 @@ class SphericalVector : public Method { FunctionSpace source_; FunctionSpace target_; -#if ATLAS_HAVE_EIGEN std::shared_ptr complexWeights_; std::shared_ptr realWeights_; -#endif + }; +#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 From fd77d76a3fc4b2dccff10e753f7afae0d665f273 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 12 Dec 2023 13:42:41 +0000 Subject: [PATCH 21/26] Disable OpenMP for older intel-classic compiler (< intel/2022.2) --- .../method/sphericalvector/SphericalVector.cc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 4f7da3481..49cbc6483 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -40,6 +40,16 @@ 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 From 3686d844643b80f3b418180459c21c8aa67b9848 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 12 Dec 2023 15:12:07 +0100 Subject: [PATCH 22/26] Enable test with CONDITION statement --- src/tests/interpolation/CMakeLists.txt | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 429323cd4..29d220590 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -122,15 +122,12 @@ ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -if ( atlas_HAVE_EIGEN ) - 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() - -endif() From 09af7e7a6805dc632fb7b96ffb4aeafd7edb6b7f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 12 Dec 2023 15:10:08 +0100 Subject: [PATCH 23/26] Revert whitespace changes --- src/atlas/CMakeLists.txt | 8 ++++---- src/atlas/interpolation/method/Method.cc | 2 +- src/atlas/interpolation/method/MethodFactory.cc | 1 + 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 1aeccbd62..3636d74d6 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 @@ -868,7 +868,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/Method.cc b/src/atlas/interpolation/method/Method.cc index 25701dbe6..fb7f90221 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -374,7 +374,7 @@ void Method::do_execute(const FieldSet& fieldsSource, FieldSet& fieldsTarget, Me ATLAS_ASSERT(N == fieldsTarget.size()); for (idx_t i = 0; i < fieldsSource.size(); ++i) { - Method::do_execute(fieldsSource[i], fieldsTarget[i], metadata); + Method::do_execute(fieldsSource[i], fieldsTarget[i], metadata); } } diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index e666e2bf2..4baa44c37 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -26,6 +26,7 @@ #include "unstructured/FiniteElement.h" #include "unstructured/UnstructuredBilinearLonLat.h" + namespace atlas { namespace interpolation { From 1212789fc107e3d42c8832a4ab55a594c9d2c70d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 12 Dec 2023 16:39:54 +0100 Subject: [PATCH 24/26] Make greatCircleCourse private before moving to eckit --- .../method/sphericalvector/SphericalVector.cc | 6 ++- src/atlas/util/Geometry.cc | 42 ++++++++++++++++++- src/atlas/util/Geometry.h | 22 ++++++++++ src/atlas/util/UnitSphere.h | 36 ---------------- src/tests/util/test_unitsphere.cc | 8 ++-- 5 files changed, 72 insertions(+), 42 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 49cbc6483..ad0311de0 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -26,7 +26,7 @@ #include "atlas/runtime/Exception.h" #include "atlas/runtime/Trace.h" #include "atlas/util/Constants.h" -#include "atlas/util/UnitSphere.h" +#include "atlas/util/Geometry.h" #include "eckit/linalg/Triplet.h" @@ -161,13 +161,15 @@ void SphericalVector::do_setup(const FunctionSpace& source, 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 = util::greatCircleCourse(sourceLonLat, targetLonLat); + const auto alpha = unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); const auto deltaAlpha = (alpha.first - alpha.second) * util::Constants::degreesToRadians(); diff --git a/src/atlas/util/Geometry.cc b/src/atlas/util/Geometry.cc index 7d8e07797..9b107b74b 100644 --- a/src/atlas/util/Geometry.cc +++ b/src/atlas/util/Geometry.cc @@ -8,12 +8,14 @@ * nor does it submit to any jurisdiction. */ +#include "atlas/util/Geometry.h" + #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 +32,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..a17f3e457 100644 --- a/src/atlas/util/Geometry.h +++ b/src/atlas/util/Geometry.h @@ -27,6 +27,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 +50,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 +79,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 +95,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 +128,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/atlas/util/UnitSphere.h b/src/atlas/util/UnitSphere.h index 1f38e455d..a3d424b77 100644 --- a/src/atlas/util/UnitSphere.h +++ b/src/atlas/util/UnitSphere.h @@ -27,42 +27,6 @@ namespace util { using eckit::geometry::UnitSphere; -/// @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 -/// -inline std::pair greatCircleCourse(const Point2& lonLat1, - const Point2& lonLat2) { - - const auto lambda1 = lonLat1[0] * Constants::degreesToRadians(); - const auto lambda2 = lonLat2[0] * Constants::degreesToRadians(); - const auto phi1 = lonLat1[1] * Constants::degreesToRadians(); - const auto phi2 = lonLat2[1] * 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 * Constants::radiansToDegrees(), - alpha2 * Constants::radiansToDegrees()); -}; - //------------------------------------------------------------------------------------------------------ } // namespace util diff --git a/src/tests/util/test_unitsphere.cc b/src/tests/util/test_unitsphere.cc index 523fcc5f9..d390334a5 100644 --- a/src/tests/util/test_unitsphere.cc +++ b/src/tests/util/test_unitsphere.cc @@ -6,7 +6,7 @@ */ #include "atlas/util/Point.h" -#include "atlas/util/UnitSphere.h" +#include "atlas/util/Geometry.h" #include "tests/AtlasTestEnvironment.h" @@ -17,6 +17,8 @@ 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.); @@ -27,8 +29,8 @@ CASE("great-circle course") { const auto targetCourse3 = 0.; const auto targetCourse4 = 180.; - const auto[ course1, course2 ] = greatCircleCourse(point1, point2); - const auto[ course3, course4 ] = greatCircleCourse(point3, point4); + 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); From 98ab54a18e9f6d485c34e08d0757a70097b31a2f Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 13 Dec 2023 16:56:30 +0100 Subject: [PATCH 25/26] Fix header includes --- src/atlas/util/Geometry.cc | 2 ++ src/atlas/util/Geometry.h | 1 + src/atlas/util/UnitSphere.h | 6 ------ 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/src/atlas/util/Geometry.cc b/src/atlas/util/Geometry.cc index 9b107b74b..074d72d34 100644 --- a/src/atlas/util/Geometry.cc +++ b/src/atlas/util/Geometry.cc @@ -10,6 +10,8 @@ #include "atlas/util/Geometry.h" +#include + #include "eckit/geometry/Point2.h" #include "eckit/geometry/Point3.h" diff --git a/src/atlas/util/Geometry.h b/src/atlas/util/Geometry.h index a17f3e457..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" diff --git a/src/atlas/util/UnitSphere.h b/src/atlas/util/UnitSphere.h index a3d424b77..c1458468f 100644 --- a/src/atlas/util/UnitSphere.h +++ b/src/atlas/util/UnitSphere.h @@ -10,12 +10,6 @@ #pragma once -#include -#include - -#include "atlas/util/Constants.h" -#include "atlas/util/Point.h" - #include "eckit/geometry/UnitSphere.h" //------------------------------------------------------------------------------------------------------ From 96b28b2682d36b3f6d5552cdf7b2787388a90e00 Mon Sep 17 00:00:00 2001 From: odlomax Date: Thu, 14 Dec 2023 17:44:30 +0000 Subject: [PATCH 26/26] Addressed reviewer comments. --- .../method/sphericalvector/SphericalVector.cc | 41 ++++++++++++++----- .../method/sphericalvector/SphericalVector.h | 15 ++----- 2 files changed, 35 insertions(+), 21 deletions(-) diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index ad0311de0..483ce8634 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -6,11 +6,15 @@ */ #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" @@ -19,17 +23,12 @@ #include "atlas/interpolation/Cache.h" #include "atlas/interpolation/Interpolation.h" #include "atlas/interpolation/method/MethodFactory.h" -#include "atlas/interpolation/method/sphericalvector/SphericalVector.h" -#include "atlas/linalg/sparse.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/linalg/Triplet.h" - namespace atlas { namespace interpolation { namespace method { @@ -189,6 +188,11 @@ void SphericalVector::do_setup(const FunctionSpace& source, 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, @@ -239,16 +243,32 @@ void SphericalVector::do_execute(const Field& sourceField, Field& targetField, 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); - } else if (sourceField.rank() == 3) { + return; + } + + if (sourceField.rank() == 3) { interpolate_vector_field(sourceField, targetField); - } else { - ATLAS_NOTIMPLEMENTED; + return; } + + ATLAS_NOTIMPLEMENTED; } template @@ -271,7 +291,9 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, matrixMultiply(sourceView, targetView, horizontalComponent, *complexWeights_); return; - } else if (sourceField.variables() == 3) { + } + + if (sourceField.variables() == 3) { const auto horizontalAndVerticalComponent = [&]( const auto& sourceVars, auto& targetVars, const auto& complexWeight, @@ -282,7 +304,6 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField, matrixMultiply(sourceView, targetView, horizontalAndVerticalComponent, *complexWeights_, *realWeights_); - return; } diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index 5807fd136..91ce2a1a2 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -19,7 +19,7 @@ #include "atlas/functionspace/FunctionSpace.h" #include "atlas/interpolation/method/Method.h" #include "atlas/linalg/sparse.h" -#include "eckit/config/Configuration.h" + namespace atlas { namespace interpolation { @@ -51,10 +51,7 @@ class SphericalVector : public Method { /// Note: This method only works with matrix-based interpolation /// schemes. /// - SphericalVector(const Config& config) : Method(config) { - const auto& conf = dynamic_cast(config); - interpolationScheme_ = conf.getSubConfiguration("scheme"); - } + SphericalVector(const Config& config); ~SphericalVector() override {} void print(std::ostream&) const override; @@ -68,13 +65,9 @@ class SphericalVector : public Method { void do_execute_adjoint(FieldSet& sourceFieldSet, const FieldSet& targetFieldSet, - Metadata& metadata) const override { - ATLAS_NOTIMPLEMENTED; - } + Metadata& metadata) const override; void do_execute_adjoint(Field& sourceField, const Field& targetField, - Metadata& metadata) const override { - ATLAS_NOTIMPLEMENTED; - } + Metadata& metadata) const override; private: template