diff --git a/src/atlas/interpolation/method/Method.cc b/src/atlas/interpolation/method/Method.cc index fb7f90221..3a2349a1d 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -301,7 +301,7 @@ void Method::setup(const FunctionSpace& source, const FunctionSpace& target) { ATLAS_TRACE("atlas::interpolation::method::Method::setup(FunctionSpace, FunctionSpace)"); this->do_setup(source, target); - if (adjoint_) { + if (adjoint_ && target.size() > 0) { Matrix tmp(*matrix_); // if interpolation is matrix free then matrix->nonZeros() will be zero. @@ -346,14 +346,14 @@ Method::Metadata Method::execute(const Field& source, Field& target) const { } Method::Metadata Method::execute_adjoint(FieldSet& source, const FieldSet& target) const { - ATLAS_TRACE("atlas::interpolation::method::Method::execute(FieldSet, FieldSet)"); + ATLAS_TRACE("atlas::interpolation::method::Method::execute_adjoint(FieldSet, FieldSet)"); Metadata metadata; this->do_execute_adjoint(source, target, metadata); return metadata; } Method::Metadata Method::execute_adjoint(Field& source, const Field& target) const { - ATLAS_TRACE("atlas::interpolation::method::Method::execute(Field, Field)"); + ATLAS_TRACE("atlas::interpolation::method::Method::execute_adjoint(Field, Field)"); Metadata metadata; this->do_execute_adjoint(source, target, metadata); return metadata; @@ -438,8 +438,8 @@ void Method::do_execute_adjoint(Field& src, const Field& tgt, Metadata&) const { throw_NotImplemented("Adjoint Interpolation does not work for fields that have missing data. ", Here()); } - if (matrix_transpose_.empty()) { - throw_AssertionFailed("Need to set 'adjoint coefficients' to true in config for adjoint interpolation to work"); + if (!adjoint_) { + throw_AssertionFailed("Need to set 'adjoint' to true in config for adjoint interpolation to work"); } if (src.datatype().kind() == array::DataType::KIND_REAL64) { diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 7e29f41c1..ab5f573d8 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -72,7 +72,7 @@ void SphericalVector::do_setup(const FunctionSpace& source, // whereas eckit does not. auto complexTriplets = ComplexTriplets(nNonZeros); auto realTriplets = RealTriplets(nNonZeros); - + const auto sourceLonLatsView = array::make_view(source_.lonlat()); const auto targetLonLatsView = array::make_view(target_.lonlat()); const auto unitSphere = geometry::UnitSphere{}; @@ -128,6 +128,11 @@ void SphericalVector::do_execute(const FieldSet& sourceFieldSet, void SphericalVector::do_execute(const Field& sourceField, Field& targetField, Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); + + if (targetField.size() == 0) { + return; + } + const auto fieldType = sourceField.metadata().getString("type", ""); if (fieldType != "vector") { auto metadata = Metadata(); @@ -162,6 +167,10 @@ void SphericalVector::do_execute_adjoint(Field& sourceField, ATLAS_TRACE( "atlas::interpolation::method::SphericalVector::do_execute_adjoint()"); + if (targetField.size() == 0) { + return; + } + const auto fieldType = sourceField.metadata().getString("type", ""); if (fieldType != "vector") { auto metadata = Metadata(); @@ -183,9 +192,6 @@ template void SphericalVector::interpolate_vector_field(const Field& sourceField, Field& targetField, const MatMul& matMul) { - if (targetField.size() == 0) { - return; - } ATLAS_ASSERT_MSG(sourceField.variables() == 2 || sourceField.variables() == 3, "Vector field can only have 2 or 3 components."); diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index b40436f97..c9376d79e 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -421,7 +421,7 @@ void StructuredInterpolation2D::setup( const FunctionSpace& source ) { } // fill sparse matrix - if( failed_points.empty() ) { + if( failed_points.empty() && out_npts_) { idx_t inp_npts = source.size(); Matrix A( out_npts_, inp_npts, triplets ); setMatrix(A); diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index 4717c5a19..95281a880 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -58,10 +58,10 @@ double vortexVertical(double lon, double lat) { void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { const auto functionSpace = fieldSet[0].functionspace(); - const auto structuredColums = functionspace::StructuredColumns(functionSpace); + const auto structuredColumns = functionspace::StructuredColumns(functionSpace); const auto nodeColumns = functionspace::NodeColumns(functionSpace); const auto mesh = - structuredColums ? Mesh(structuredColums.grid()) : nodeColumns.mesh(); + structuredColumns ? Mesh(structuredColumns.grid()) : nodeColumns.mesh(); const auto gmshConfig = Config("coordinates", "xyz") | Config("ghost", true) | Config("info", true); @@ -79,10 +79,16 @@ const auto generateNodeColums(const std::string& gridName, return functionspace::NodeColumns(mesh); } +// Helper function to create part-empty PointCloud +const auto generateEmptyPointCloud() { + const auto functionSpace = functionspace::PointCloud(std::vector{}); + return functionSpace; +} + // Helper struct to key different Functionspaces to strings struct FunctionSpaceFixtures { static const FunctionSpace& get(const std::string& fixture) { - static const auto functionSpaces = + static auto functionSpaces = std::map{ {"cubedsphere_mesh", generateNodeColums("CS-LFR-48", "cubedsphere_dual")}, @@ -92,7 +98,8 @@ struct FunctionSpaceFixtures { {"structured_columns_lowres", functionspace::StructuredColumns(Grid("O24"), option::halo(1))}, {"structured_columns_hires", - functionspace::StructuredColumns(Grid("O96"), option::halo(1))}}; + functionspace::StructuredColumns(Grid("O96"), option::halo(1))}, + {"empty_point_cloud", generateEmptyPointCloud()}}; return functionSpaces.at(fixture); } }; @@ -109,14 +116,16 @@ struct FieldSpecFixtures { } }; -// Helper stcut to key different interpolation schemes to strings +// Helper struct to key different interpolation schemes to strings struct InterpSchemeFixtures { static const Config& get(const std::string& fixture) { static const auto cubedsphereBilinear = - option::type("cubedsphere-bilinear"); - static const auto finiteElement = option::type("finite-element"); - static const auto structuredLinear = - option::type("structured-linear2D") | option::halo(1); + option::type("cubedsphere-bilinear") | Config("adjoint", true); + static const auto finiteElement = + option::type("finite-element") | Config("adjoint", true); + static const auto structuredLinear = option::type("structured-linear2D") | + option::halo(1) | + Config("adjoint", true); static const auto sphericalVector = option::type("spherical-vector") | Config("adjoint", true); @@ -218,40 +227,45 @@ void testInterpolation(const Config& config) { targetFunctionSpace.createField(errorFieldSpec))); errorView.assign(0.); - auto maxError = 0.; - ArrayForEach<0>::apply( - std::tie(targetLonLat, targetView, errorView), - [&](auto&& lonLat, auto&& targetColumn, auto&& errorColumn) { - const auto calcError = [&](auto&& targetElem, auto&& errorElem) { - auto trueValue = std::vector(targetElem.size()); - std::tie(trueValue[0], trueValue[1]) = - vortexHorizontal(lonLat(0), lonLat(1)); - if (targetElem.size() == 3) { - trueValue[2] = vortexVertical(lonLat(0), lonLat(1)); - } - - auto errorSqrd = 0.; - for (auto k = 0; k < targetElem.size(); ++k) { - errorSqrd += - (targetElem(k) - trueValue[k]) * (targetElem(k) - trueValue[k]); + if (config.has("tol")) { + auto maxError = 0.; + ArrayForEach<0>::apply( + std::tie(targetLonLat, targetView, errorView), + [&](auto&& lonLat, auto&& targetColumn, auto&& errorColumn) { + const auto calcError = [&](auto&& targetElem, auto&& errorElem) { + auto trueValue = std::vector(targetElem.size()); + std::tie(trueValue[0], trueValue[1]) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (targetElem.size() == 3) { + trueValue[2] = vortexVertical(lonLat(0), lonLat(1)); + } + + auto errorSqrd = 0.; + for (auto k = 0; k < targetElem.size(); ++k) { + errorSqrd += (targetElem(k) - trueValue[k]) * + (targetElem(k) - trueValue[k]); + } + + errorElem = std::sqrt(errorSqrd); + maxError = std::max(maxError, static_cast(errorElem)); + }; + + if constexpr (Rank == 2) { + calcError(targetColumn, errorColumn); } + else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), + calcError); + } + }); - 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")); + EXPECT_APPROX_EQ(maxError, 0., config.getDouble("tol")); + } - gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); - gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); + if (config.has("file_id")) { + gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); + gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); + } // Adjoint test auto targetAdjoint = targetFunctionSpace.createField(fieldSpec); @@ -275,8 +289,11 @@ void testInterpolation(const Config& config) { constexpr auto tinyNum = 1e-13; const auto targetDotTarget = dotProduct(targetView, targetView); const auto sourceDotSourceAdjoint = dotProduct(sourceView, sourceAdjointView); - const auto dotProdRatio = targetDotTarget / sourceDotSourceAdjoint; - EXPECT_APPROX_EQ(dotProdRatio, 1., tinyNum); + + if (targetFunctionSpace.size() > 0) { + const auto dotProdRatio = targetDotTarget / sourceDotSourceAdjoint; + EXPECT_APPROX_EQ(dotProdRatio, 1., tinyNum); + } } CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { @@ -347,6 +364,44 @@ CASE("structured columns vector interpolation (2d-field, 2-vector, hi-res)") { testInterpolation((config)); } +CASE("cubed sphere (spherical vector) to empty point cloud") { + const auto config = + Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "empty_point_cloud") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "cubedsphere_bilinear_spherical"); + + testInterpolation((config)); +} + +CASE("cubed sphere to empty point cloud") { + const auto config = + Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "empty_point_cloud") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "cubedsphere_bilinear"); + + testInterpolation((config)); +} + +CASE("structured columns to empty point cloud") { + const auto config = Config("source_fixture", "structured_columns") + .set("target_fixture", "empty_point_cloud") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear"); + + testInterpolation((config)); +} + +CASE("finite element to empty point cloud") { + const auto config = Config("source_fixture", "gaussian_mesh") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "finite_element"); + + testInterpolation((config)); +} + } // namespace test } // namespace atlas