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

Allow zero-sized interpolation target functionspace #206

Merged
merged 4 commits into from
Jun 7, 2024
Merged

Allow zero-sized interpolation target functionspace #206

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f5da33f
Added failing tests.
odlomax Jun 6, 2024
6c91c4a
Fixed StructuredInterpolation2D.
odlomax Jun 6, 2024
4f63ada
Fixed SphericalVector.
odlomax Jun 6, 2024
5fada06
Removed redundant headers.
odlomax Jun 6, 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
10 changes: 5 additions & 5 deletions src/atlas/interpolation/method/Method.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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) {
Expand Down
14 changes: 10 additions & 4 deletions src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double, 2>(source_.lonlat());
const auto targetLonLatsView = array::make_view<double, 2>(target_.lonlat());
const auto unitSphere = geometry::UnitSphere{};
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand All @@ -183,9 +192,6 @@ template <typename MatMul>
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.");
Expand Down
2 changes: 1 addition & 1 deletion src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -421,7 +421,7 @@ void StructuredInterpolation2D<Kernel>::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);
Expand Down
139 changes: 97 additions & 42 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 @@ -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);
Expand All @@ -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<PointXY>{});
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<std::string_view, FunctionSpace>{
{"cubedsphere_mesh",
generateNodeColums("CS-LFR-48", "cubedsphere_dual")},
Expand All @@ -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);
}
};
Expand All @@ -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);
Expand Down Expand Up @@ -218,40 +227,45 @@ void testInterpolation(const Config& config) {
targetFunctionSpace.createField<double>(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<double>(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<double>(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<double>(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<double>(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<double>(fieldSpec);
Expand All @@ -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)") {
Expand Down Expand Up @@ -347,6 +364,44 @@ CASE("structured columns vector interpolation (2d-field, 2-vector, hi-res)") {
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");

testInterpolation<Rank2dField>((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<Rank2dField>((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<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));
}

} // namespace test
} // namespace atlas

Expand Down
Loading