Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bug fixes to vector interpolation with StructuredColumns and spherical vector interpolation #222

Merged
merged 2 commits into from
Sep 16, 2024
Merged

Bug fixes to vector interpolation with StructuredColumns and spherical vector interpolation #222

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9268d27
bug fixes
MarekWlasak Sep 12, 2024
f8e805e
removing commented out code
MarekWlasak Sep 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions src/atlas/functionspace/detail/StructuredColumns.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}
Expand Down Expand Up @@ -726,7 +726,7 @@ struct FixupHaloForVectors {
template <typename DATATYPE>
void apply(Field& field) {
std::string type = field.metadata().getString("type", "scalar");
if (type == "vector ") {
if (type == "vector") {
ATLAS_NOTIMPLEMENTED;
}
}
Expand Down Expand Up @@ -770,12 +770,16 @@ struct FixupHaloForVectors<3> {
template <typename DATATYPE>
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<DATATYPE, RANK>(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);
}
Expand All @@ -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);
}
Expand Down
2 changes: 1 addition & 1 deletion src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ class ComplexMatrixMultiply {
// tinyNum ~= 2.3e-13 for double.
constexpr auto tinyNum = 1024 * std::numeric_limits<Real>::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 = [&]() {
Expand Down
151 changes: 107 additions & 44 deletions src/tests/interpolation/test_interpolation_spherical_vector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -126,20 +134,25 @@ 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);

static const auto interpSchemes = std::map<std::string_view, Config>{
{"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);
}
};
Expand Down Expand Up @@ -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"));
Expand Down Expand Up @@ -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")
Expand All @@ -308,7 +322,7 @@ CASE("cubed sphere vector interpolation (3d-field, 2-vector)") {
testInterpolation<Rank3dField>((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")
Expand All @@ -320,6 +334,36 @@ CASE("cubed sphere vector interpolation (3d-field, 3-vector)") {
testInterpolation<Rank3dField>((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<Rank2dField>((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<Rank2dField>((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<Rank2dField>((config));
}

CASE("finite element vector interpolation (2d-field, 2-vector)") {
const auto config = Config("source_fixture", "gaussian_mesh")
.set("target_fixture", "cubedsphere_mesh")
Expand All @@ -331,60 +375,77 @@ CASE("finite element vector interpolation (2d-field, 2-vector)") {
testInterpolation<Rank2dField>((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<Rank2dField>((config));
testInterpolation<Rank3dField>((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<Rank2dField>((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<Rank3dField>((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<Rank2dField>((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<Rank2dField>((config));
testInterpolation<Rank3dField>((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<Rank2dField>((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")
Expand All @@ -393,11 +454,13 @@ CASE("structured columns to empty point cloud") {
testInterpolation<Rank2dField>((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<Rank2dField>((config));
}
Expand Down
Loading