diff --git a/src/atlas/functionspace/detail/StructuredColumns.cc b/src/atlas/functionspace/detail/StructuredColumns.cc index 653ea10ab..54b215feb 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.cc +++ b/src/atlas/functionspace/detail/StructuredColumns.cc @@ -540,7 +540,7 @@ const util::PartitionPolygon& StructuredColumns::polygon(idx_t halo) const { } if (not polygons_[halo]) { if (halo > this->halo()) { - throw_Exception("StructuredColumsn does not contain a halo of size " + std::to_string(halo) + ".", Here()); + throw_Exception("StructuredColumns does not contain a halo of size " + std::to_string(halo) + ".", Here()); } polygons_[halo].reset(new grid::StructuredPartitionPolygon(*this, halo)); } @@ -726,7 +726,7 @@ struct FixupHaloForVectors { template void apply(Field& field) { std::string type = field.metadata().getString("type", "scalar"); - if (type == "vector ") { + if (type == "vector") { ATLAS_NOTIMPLEMENTED; } } @@ -770,12 +770,16 @@ struct FixupHaloForVectors<3> { template void apply(Field& field) { std::string type = field.metadata().getString("type", "scalar"); + + // if levels not setup in functionspace use levels from field + idx_t k_end = (fs.k_begin() == 0) && (fs.k_end() == 0) ? field.levels() : fs.k_end(); + if (type == "vector") { auto array = array::make_view(field); for (idx_t j = fs.j_begin_halo(); j < 0; ++j) { for (idx_t i = fs.i_begin_halo(j); i < fs.i_end_halo(j); ++i) { idx_t n = fs.index(i, j); - for (idx_t k = fs.k_begin(); k < fs.k_end(); ++k) { + for (idx_t k = fs.k_begin(); k < k_end; ++k) { array(n, k, XX) = -array(n, k, XX); array(n, k, YY) = -array(n, k, YY); } @@ -784,7 +788,7 @@ struct FixupHaloForVectors<3> { for (idx_t j = fs.grid().ny(); j < fs.j_end_halo(); ++j) { for (idx_t i = fs.i_begin_halo(j); i < fs.i_end_halo(j); ++i) { idx_t n = fs.index(i, j); - for (idx_t k = fs.k_begin(); k < fs.k_end(); ++k) { + for (idx_t k = fs.k_begin(); k < k_end; ++k) { array(n, k, XX) = -array(n, k, XX); array(n, k, YY) = -array(n, k, YY); } diff --git a/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h b/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h index 34cfdacb2..8d15d035a 100644 --- a/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h +++ b/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h @@ -70,7 +70,7 @@ class ComplexMatrixMultiply { // tinyNum ~= 2.3e-13 for double. constexpr auto tinyNum = 1024 * std::numeric_limits::epsilon(); const auto complexMagnitude = std::abs(complexRowIter.value()); - const auto realValue = realRowIter.value(); + const auto realValue = std::abs(realRowIter.value()); const auto error = std::abs(complexMagnitude - realValue); const auto printError = [&]() { diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc index 95281a880..16dbdabf9 100644 --- a/src/tests/interpolation/test_interpolation_spherical_vector.cc +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -95,6 +95,14 @@ struct FunctionSpaceFixtures { {"gaussian_mesh", generateNodeColums("O48", "structured")}, {"structured_columns", functionspace::StructuredColumns(Grid("O48"), option::halo(1))}, + {"structured_columns_classic", + functionspace::StructuredColumns(Grid("F48"), option::halo(1))}, + {"structured_columns_classic_halo2", + functionspace::StructuredColumns(Grid("F48"), option::halo(2))}, + {"structured_columns_classic_highres_halo2", + functionspace::StructuredColumns(Grid("F96"), option::halo(2))}, + {"structured_columns_halo2", + functionspace::StructuredColumns(Grid("O48"), option::halo(2))}, {"structured_columns_lowres", functionspace::StructuredColumns(Grid("O24"), option::halo(1))}, {"structured_columns_hires", @@ -126,7 +134,9 @@ struct InterpSchemeFixtures { static const auto structuredLinear = option::type("structured-linear2D") | option::halo(1) | Config("adjoint", true); - + static const auto structuredCubic = option::type("structured-bicubic") | + option::halo(2) | + Config("adjoint", true); static const auto sphericalVector = option::type("spherical-vector") | Config("adjoint", true); @@ -134,12 +144,15 @@ struct InterpSchemeFixtures { {"cubedsphere_bilinear", cubedsphereBilinear}, {"finite_element", finiteElement}, {"structured_linear", structuredLinear}, + {"structured_cubic", structuredCubic}, {"cubedsphere_bilinear_spherical", sphericalVector | Config("scheme", cubedsphereBilinear)}, {"finite_element_spherical", sphericalVector | Config("scheme", finiteElement)}, {"structured_linear_spherical", - sphericalVector | Config("scheme", structuredLinear)}}; + sphericalVector | Config("scheme", structuredLinear)}, + {"structured_cubic_spherical", + sphericalVector | Config("scheme", structuredCubic)}}; return interpSchemes.at(fixture); } }; @@ -254,9 +267,9 @@ void testInterpolation(const Config& config) { calcError(targetColumn, errorColumn); } else if constexpr (Rank == 3) { - ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), - calcError); - } + ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), + calcError); + } }); EXPECT_APPROX_EQ(maxError, 0., config.getDouble("tol")); @@ -296,7 +309,8 @@ void testInterpolation(const Config& config) { } } -CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { + +CASE("cubed sphere CS-LFR-48 vector interpolation (3d-field, 2-vector)") { const auto config = Config("source_fixture", "cubedsphere_mesh") .set("target_fixture", "gaussian_mesh") @@ -308,7 +322,7 @@ CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { testInterpolation((config)); } -CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { +CASE("cubed sphere CS-LFR-48 vector interpolation (3d-field, 3-vector)") { const auto config = Config("source_fixture", "cubedsphere_mesh") .set("target_fixture", "gaussian_mesh") @@ -320,6 +334,36 @@ CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { testInterpolation((config)); } +CASE("cubed sphere CS-LFR-48 (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 CS-LFR-48 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("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)); +} + CASE("finite element vector interpolation (2d-field, 2-vector)") { const auto config = Config("source_fixture", "gaussian_mesh") .set("target_fixture", "cubedsphere_mesh") @@ -331,60 +375,77 @@ CASE("finite element vector interpolation (2d-field, 2-vector)") { testInterpolation((config)); } -CASE("structured columns vector interpolation (2d-field, 2-vector)") { - const auto config = Config("source_fixture", "structured_columns") - .set("target_fixture", "cubedsphere_mesh") - .set("field_spec_fixture", "2vector") - .set("interp_fixture", "structured_linear_spherical") - .set("file_id", "spherical_vector_sc") - .set("tol", 0.00017); +CASE("structured columns F48 cubic vector spherical interpolation (3d-field, 2-vector)") { + const auto config = + Config("source_fixture", "structured_columns_classic_halo2") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_cubic_spherical") + .set("file_id", "spherical_cubic_vector_classic_sc") + .set("tol", 0.0000082); - testInterpolation((config)); + testInterpolation((config)); } -CASE("structured columns vector interpolation (2d-field, 2-vector, low-res)") { - const auto config = Config("source_fixture", "structured_columns_lowres") - .set("target_fixture", "gaussian_mesh") - .set("field_spec_fixture", "2vector") - .set("interp_fixture", "structured_linear_spherical") - .set("file_id", "spherical_vector_sc_lr") - .set("tol", 0.00056); +CASE("structured columns F96 cubic vector spherical interpolation (2d-field, 2-vector)") { + const auto config = + Config("source_fixture", "structured_columns_classic_highres_halo2") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_cubic_spherical") + .set("file_id", "spherical_2D_cubic_vector_highres_classic_sc") + .set("tol", 0.000000425); testInterpolation((config)); } -CASE("structured columns vector interpolation (2d-field, 2-vector, hi-res)") { - const auto config = Config("source_fixture", "structured_columns_hires") +CASE("structured columns F96 cubic vector spherical interpolation (3d-field, 2-vector)") { + const auto config = + Config("source_fixture", "structured_columns_classic_highres_halo2") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_cubic_spherical") + .set("file_id", "spherical_3D_cubic_vector_highres_classic_sc") + .set("tol", 0.00000085); + + testInterpolation((config)); +} + +CASE("structured columns O24 linear vector interpolation (2d-field, 2-vector)") { + const auto config = Config("source_fixture", "structured_columns_lowres") .set("target_fixture", "gaussian_mesh") .set("field_spec_fixture", "2vector") .set("interp_fixture", "structured_linear_spherical") - .set("file_id", "spherical_vector_sc_hr") - .set("tol", 0.000044); + .set("file_id", "spherical_vector_sc_lr") + .set("tol", 0.00056); 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"); +CASE("structured columns O48 cubic vector spherical interpolation (3d-field, 2-vector)") { + const auto config = + Config("source_fixture", "structured_columns_halo2") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_cubic_spherical") + .set("file_id", "spherical_cubic_vector_sc3") + .set("tol", 0.000007); - testInterpolation((config)); + 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"); +CASE("structured columns O48 linear vector interpolation (2d-field, 2-vector)") { + const auto config = Config("source_fixture", "structured_columns") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc") + .set("tol", 0.00017); testInterpolation((config)); } -CASE("structured columns to empty point cloud") { +CASE("structured columns O48 to empty point cloud") { const auto config = Config("source_fixture", "structured_columns") .set("target_fixture", "empty_point_cloud") .set("field_spec_fixture", "2vector") @@ -393,11 +454,13 @@ CASE("structured columns to empty point cloud") { testInterpolation((config)); } -CASE("finite element to empty point cloud") { - const auto config = Config("source_fixture", "gaussian_mesh") - .set("target_fixture", "cubedsphere_mesh") +CASE("structured columns O96 vector interpolation (2d-field, 2-vector, hi-res)") { + const auto config = Config("source_fixture", "structured_columns_hires") + .set("target_fixture", "gaussian_mesh") .set("field_spec_fixture", "2vector") - .set("interp_fixture", "finite_element"); + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_hr") + .set("tol", 0.000044); testInterpolation((config)); }