From 346f7f1c904abfb0bd38590c249a565ed4177563 Mon Sep 17 00:00:00 2001 From: odlomax Date: Tue, 5 Dec 2023 10:35:14 +0000 Subject: [PATCH] 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)); +} } }